RHESSysIOinR_Coweeta_WS32_74 Multi Site Validation Field Data_Step1
Carlos Quintero
6/12/2024
Coweeta Hydrologic Laboratory Watershed 32 using Static, SSURGO, and RSS soil inputs.
The following R Markdown script prepares and runs the Regional Hydro-Ecologic Simulation System (RHESSYS) within the R environment. Modified from RHESSysIOinR example scripts created by Will Burke and Ryan Bart. Includes CLHS code snippets and SDA code from Dylan Beaudette.
Windows Subsystem via Linux must be installed before RHESSys will run in a Windows environment. GRASS8 must be installed to run the spatial preprocessing portion of the script.
Code has been modified to run RHESSys 7.4 and 5.18. 7.4 vegetation processing differs slightly from 5.18. Various template differences exist between the two versions of RHESSys and template files must be carefully reviewed if moving between versions. This script is currently set up to run RHESSys 7.4
knitr::opts_chunk$set(warning = FALSE, message = FALSE)
#install.packages("devtools")
#devtools::install_github("ropensci/FedData")
library(aqp)## This is aqp 1.41library(daymetr)
library(clhs)## Registered S3 methods overwritten by 'tibble':
## method from
## format.tbl pillar
## print.tbl pillarlibrary(EcoHydRology)## Loading required package: operators##
## Attaching package: 'operators'## The following objects are masked from 'package:base':
##
## options, strrep## Loading required package: topmodel## Loading required package: DEoptim## Loading required package: parallel##
## DEoptim package
## Differential Evolution algorithm in R
## Authors: D. Ardia, K. Mullen, B. Peterson and J. Ulrich## Loading required package: XMLlibrary(elevatr)
library(ezknitr)
library(FedData)##
## Attaching package: 'FedData'## The following object is masked from 'package:operators':
##
## %>%library(foreach)
library(foreign)
library(ggplot2)
library(ggpubr)##
## Attaching package: 'ggpubr'
##
## The following object is masked from 'package:operators':
##
## %>%library(hydroGOF)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric##
## Attaching package: 'hydroGOF'## The following object is masked from 'package:topmodel':
##
## NSefflibrary(imputeTS)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo##
## Attaching package: 'imputeTS'## The following object is masked from 'package:zoo':
##
## na.locf## The following object is masked from 'package:operators':
##
## %>%library(knitr)
library(latticeExtra)## Loading required package: lattice##
## Attaching package: 'latticeExtra'## The following object is masked from 'package:ggplot2':
##
## layerlibrary(lubridate)##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionlibrary(Metrics)##
## Attaching package: 'Metrics'## The following objects are masked from 'package:hydroGOF':
##
## mae, mse, rmselibrary(raster)## Loading required package: sp##
## Attaching package: 'raster'## The following objects are masked from 'package:aqp':
##
## metadata, metadata<-library(rasterVis)
library(readxl)
library(reshape2)
library(rgdal)## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
##
## rgdal: version: 1.5-32, (SVN revision 1176)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.1, released 2021/12/27
## Path to GDAL shared files: C:/Users/Carlos/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/Carlos/Documents/R/win-library/4.0/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-6
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.library(rgrass)## GRASS GIS interface loaded with GRASS version: GRASS 8.2.1 (2023)
## and location: Coweeta1library(rmarkdown)
library(Rmisc)## Loading required package: plyr##
## Attaching package: 'plyr'## The following object is masked from 'package:ggpubr':
##
## mutatelibrary(RHESSysPreprocessing)
library(RHESSysIOinR)##
## Attaching package: 'RHESSysIOinR'## The following object is masked from 'package:RHESSysPreprocessing':
##
## read_worldlibrary(sf)## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE##
## Attaching package: 'sf'## The following object is masked from 'package:operators':
##
## %>%library(soilDB)
library(sp)
library(stats)
library(stringi)
library(tactile)
library(terra)## terra 1.5.21##
## Attaching package: 'terra'## The following object is masked from 'package:rgdal':
##
## project## The following object is masked from 'package:knitr':
##
## spin## The following object is masked from 'package:zoo':
##
## time<-## The following object is masked from 'package:ggpubr':
##
## rotate## The following object is masked from 'package:ggplot2':
##
## arrowlibrary(testthat)##
## Attaching package: 'testthat'## The following object is masked from 'package:terra':
##
## describe## The following object is masked from 'package:operators':
##
## %>%library(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v tibble 3.0.3 v dplyr 1.0.5
## v tidyr 1.1.3 v stringr 1.4.0
## v readr 2.1.2 v forcats 0.5.1
## v purrr 0.3.4## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x forcats::%>%() masks stringr::%>%(), dplyr::%>%(), purrr::%>%(), tidyr::%>%(), tibble::%>%(), testthat::%>%(), sf::%>%(), imputeTS::%>%(), ggpubr::%>%(), FedData::%>%(), operators::%>%()
## x purrr::accumulate() masks foreach::accumulate()
## x dplyr::arrange() masks plyr::arrange()
## x terra::arrow() masks ggplot2::arrow()
## x lubridate::as.difftime() masks base::as.difftime()
## x dplyr::combine() masks aqp::combine()
## x purrr::compact() masks plyr::compact()
## x dplyr::count() masks plyr::count()
## x lubridate::date() masks base::date()
## x readr::edition_get() masks testthat::edition_get()
## x tidyr::extract() masks terra::extract(), raster::extract()
## x dplyr::failwith() masks plyr::failwith()
## x dplyr::filter() masks stats::filter()
## x dplyr::id() masks plyr::id()
## x terra::intersect() masks raster::intersect(), lubridate::intersect(), base::intersect()
## x purrr::is_null() masks testthat::is_null()
## x dplyr::lag() masks stats::lag()
## x latticeExtra::layer() masks ggplot2::layer()
## x readr::local_edition() masks testthat::local_edition()
## x dplyr::matches() masks tidyr::matches(), testthat::matches()
## x dplyr::mutate() masks plyr::mutate(), ggpubr::mutate()
## x dplyr::rename() masks plyr::rename()
## x dplyr::select() masks raster::select()
## x lubridate::setdiff() masks base::setdiff()
## x dplyr::slice() masks aqp::slice()
## x dplyr::src() masks terra::src()
## x dplyr::summarise() masks plyr::summarise()
## x dplyr::summarize() masks plyr::summarize()
## x terra::union() masks raster::union(), lubridate::union(), base::union()
## x purrr::when() masks foreach::when()library(viridisLite)
library(dplyr)
#knitr::opts_knit$set(root.dir = system.file("extdata/", package = "RHESSysIOinR"))
setwd(system.file("extdata/", package = "RHESSysIOinR"))
### Downloads most recent version of RHESSYs
# gert::git_clone(url = "https://github.com/RHESSys/RHESSys", path = "/modelfiles/rhessys", branch = "trunk")
# compile_rhessys(location = "modelfiles/")
rh_ver = dir(path = "modelfiles/rhessys/", pattern = "^rhessys\\d+",recursive = F)
getwd()## [1] "C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata"rh_path = file.path("modelfiles/rhessys", rh_ver)
## Load in necessary files
filedir <- ("C:/Users/Carlos/Documents/GitHub/RHESSys74-R-Markdown/")
watershedshapefile <- paste0(filedir, "input_files/shapefiles/watersheds/Coweeta_Hydrologic_Laboratory.shp")
streamshapefile<- paste0(filedir,"input_files/shapefiles/streams/StreamsCaldwellClipped.shp")
weirshapefile<- paste0(filedir,"input_files/shapefiles/weirs/Subbasinoutlets.shp")
roadshapefile <- paste0(filedir, "input_files/shapefiles/roads/UpdatedRoads4Reprojected32617.shp")
CWRG12 <- paste0(filedir,"input_files/precipitation/RDS-2017-0031/Data/RG12_daily_1942_2021.csv")
CWRG20 <-paste0(filedir,"input_files/precipitation/RDS-2017-0031/Data/RG20_daily_1962_2021.csv")
CS01 <- ("C:/Users/Carlos/Desktop/ORISE/Climate Data/CW/Air Temperature/CS01 OISHI/cs01_hourly_Tair_2011_2023.csv")
CS01LT<- ("C:/Users/Carlos/Desktop/ORISE/Climate Data/CW/Air Temperature/RDS-2015-0049/Data/cs01_daily_airtemp.csv")
CS21<- paste0(filedir,"input_files/temperature/RDS-2015-0042/Data/cs21_daily.csv")
CS28<- paste0(filedir,"input_files/temperature/RDS-2015-0042/Data/cs28_daily.csv")
newsmlocations<- paste0(filedir,"input_files/shapefiles/soil_moisture/CoweetaNRCSsmKML.kml")Create Local Base File if one is not available
setwd(system.file("extdata/", package = "RHESSysIOinR"))
localbasefiletext<- "101 base_station_id
275535 x_coordinate
3879165 y_coordinate
1064 z_coordinate
3.5 effective_lai
5.0 screen_height
clim/cwt_annual annual_climate_prefix
0 number_non_critical_annual_sequences
clim/cwt_monthly monthly_climate_prefix
0 number_non_critical_monthly_sequences
clim/cwtws32local daily_climate_prefix
0 number_non_critical_daily_sequences
clim/cwt_hourly hourly_climate_prefix
0 number_non_critical_hourly_sequences
clim/cwt trigger_of_dated_input
0 num_non_critical_sequences
"
write(localbasefiletext,file = "clim/cwtws32local.base")Download or find Local Climate for Watershed of Interest. Precipitation and Temperature
https://www.fs.usda.gov/rds/archive/Catalog/RDS-2017-0031
RGO6-WS32 & RG12-WS32 - Using RG12 RG20-WS7 Precipitation is in inches
CS21-WS7 CS28-WS32 Temperature is in celsius
RHESSys can use Rain (Meters),Tmin (C),Tmax (C),Tavg (C),RHmin (0-1),RHmax (0-1),RHavg (0-1),Wind (Meters/Sec),Wind Dir (Degrees)
only using Rain, Tmin and Tmax for now
“One of the most intense arctic outbreaks of the 20th century occurred January 18-22, 1985. Extremely cold temperatures affected every state east of the Rockies with three new state record lows established: -34° at Mt. Mitchell, North Carolina, -19° at Caesar’s Head, South Carolina; and -30° at Mountain Lake, Virginia.” Taken from weather.gov
lowest temp -34, all TMIN values less than -34 will be removed as erroneous values
setwd(system.file("extdata/", package = "RHESSysIOinR"))
cwtrg12precip<- read.table(file=CWRG12, header = TRUE, sep = ",")
cwtrg20precip<- read.table(file=CWRG20, header = TRUE, sep = ",")
cwtrg12precip$date <- as.Date(paste(cwtrg12precip$YEAR,cwtrg12precip$MONTH,cwtrg12precip$DAY, sep = "-"), "%Y-%m-%d")
cwtrg20precip$date <- as.Date(paste(cwtrg20precip$YEAR,cwtrg20precip$MONTH,cwtrg20precip$DAY, sep = "-"), "%Y-%m-%d")
plot(cwtrg12precip$date,cwtrg12precip$RRG12, type = "h", main = "Coweeta Rain Gauge 12 Precipitation", col = "BLUE")
plot(cwtrg20precip$date,cwtrg20precip$RRG20, type = "h", main = "Coweeta Rain Gauge 20 Precipitation", col = 4)
cwtcs21<- read.table(file=CS21, header = TRUE, sep = ",")
cwtcs28<- read.table(file=CS28, header = TRUE, sep = ",")
cwtcs21$date<- as.Date(paste(cwtcs21$YEAR,cwtcs21$MONTH,cwtcs21$DAY, sep = "-"), "%Y-%m-%d")
cwtcs28$date<- as.Date(paste(cwtcs28$YEAR,cwtcs28$MONTH,cwtcs28$DAY, sep = "-"), "%Y-%m-%d")
cwtcs21 <- as_tibble(cwtcs21)
cwtcs28 <- as_tibble(cwtcs28)
##filtering low temp in climate station files
cwtcs21<- cwtcs21 %>% mutate(TMIN= ifelse(TMIN < -34, NA,TMIN))
cwtcs28<- cwtcs28 %>% mutate(TMIN= ifelse(TMIN < -34, NA,TMIN))
#Combine Rain Gauge and Climate Station datasets
cwtrg12precip<- as_tibble(cwtrg12precip)
cwtrg20precip<- as_tibble(cwtrg20precip)
###make sure to create full date time series before left join
full_dates<- tibble(date = seq.Date(from = min(cwtrg12precip$date), to = max(cwtrg12precip$date), by = "day"))
full_dates## # A tibble: 29,066 x 1
## date
## <date>
## 1 1942-06-04
## 2 1942-06-05
## 3 1942-06-06
## 4 1942-06-07
## 5 1942-06-08
## 6 1942-06-09
## 7 1942-06-10
## 8 1942-06-11
## 9 1942-06-12
## 10 1942-06-13
## # ... with 29,056 more rowscombinedcwtclimate<- full_dates %>% left_join(cwtrg12precip, by = c("date"))
cwtrg20_selected <- cwtrg20precip %>% select(date, RRG20)
combinedcwtclimate <- combinedcwtclimate %>% left_join(cwtrg20_selected, by = c("date"))
#convert precip from inces to meters
combinedcwtclimate$RRG12 <- combinedcwtclimate$RRG12* 0.0254
combinedcwtclimate$RRG20 <- combinedcwtclimate$RRG20* 0.0254
##precip converted from inches to meters
## RH converted to 0-1 range
cwtcs21_selected <- cwtcs21 %>% mutate(RHMIN = RHMIN/100, RHMAX = RHMAX/100, RHAVG = RHAVG/100) %>% select(date, TMINcs21=TMIN, TMAXcs21=TMAX, TAVGcs21=TAVG, RHMINcs21=RHMIN, RHMAXcs21=RHMAX, RHcs21=RHAVG, WINDcs21=WIND, WDIRcs21=W_DIR)
cwtcs28_selected <- cwtcs28 %>% mutate(RHMIN = RHMIN/100, RHMAX = RHMAX/100, RHAVG = RHAVG/100) %>% select(date, TMINcs28=TMIN, TMAXcs28=TMAX, TAVGcs28=TAVG, RHMINcs28=RHMIN, RHMAXcs28=RHMAX, RHcs28=RHAVG, WINDcs28=WIND, WDIRcs28=W_DIR)
combinedcwtclimate <- combinedcwtclimate %>% left_join(cwtcs21_selected, by = c("date"))
combinedcwtclimate <- combinedcwtclimate %>% left_join(cwtcs28_selected, by = c("date"))
combinedcwtclimate## # A tibble: 29,066 x 22
## date YEAR MONTH DAY RRG12 RRG20 TMINcs21 TMAXcs21 TAVGcs21
## <date> <int> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1942-06-04 1942 6 4 1.27e-3 NA NA NA NA
## 2 1942-06-05 1942 6 5 2.54e-4 NA NA NA NA
## 3 1942-06-06 NA NA NA NA NA NA NA NA
## 4 1942-06-07 1942 6 7 9.40e-3 NA NA NA NA
## 5 1942-06-08 1942 6 8 9.40e-3 NA NA NA NA
## 6 1942-06-09 1942 6 9 1.09e-2 NA NA NA NA
## 7 1942-06-10 1942 6 10 4.19e-2 NA NA NA NA
## 8 1942-06-11 1942 6 11 8.89e-3 NA NA NA NA
## 9 1942-06-12 1942 6 12 6.60e-3 NA NA NA NA
## 10 1942-06-13 1942 6 13 5.08e-4 NA NA NA NA
## # ... with 29,056 more rows, and 13 more variables: RHMINcs21 <dbl>,
## # RHMAXcs21 <dbl>, RHcs21 <dbl>, WINDcs21 <dbl>, WDIRcs21 <dbl>,
## # TMINcs28 <dbl>, TMAXcs28 <dbl>, TAVGcs28 <dbl>, RHMINcs28 <dbl>,
## # RHMAXcs28 <dbl>, RHcs28 <dbl>, WINDcs28 <dbl>, WDIRcs28 <dbl>first_valid_index<- combinedcwtclimate %>% dplyr::mutate(row_num = row_number()) %>% filter(!is.na(RRG12)& !is.na(RRG20)& !is.na(TMINcs21)& !is.na(TMINcs28) ) %>% slice(1) %>% pull(row_num)
last_valid_index<- combinedcwtclimate %>% dplyr::mutate(row_num = row_number()) %>% filter(!is.na(RRG12)& !is.na(RRG20)& !is.na(TMINcs21)& !is.na(TMINcs28) ) %>% slice(n()) %>% pull(row_num)
first_valid_index## [1] 17382trimmed_cwtclimate <- combinedcwtclimate %>% slice(first_valid_index:last_valid_index)
trimmed_cwtclimate## # A tibble: 10,954 x 22
## date YEAR MONTH DAY RRG12 RRG20 TMINcs21 TMAXcs21 TAVGcs21
## <date> <int> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1990-01-04 1990 1 4 0.0429 0.0325 6.9 11.7 9.71
## 2 1990-01-05 1990 1 5 0.00533 0.00381 5.5 14.3 8.73
## 3 1990-01-06 1990 1 6 0.00686 0.00737 5.5 11.1 7.00
## 4 1990-01-07 1990 1 7 0.0112 0.00991 5 6.8 5.9
## 5 1990-01-08 1990 1 8 0.00584 0.00406 -0.2 13.2 6.46
## 6 1990-01-09 NA NA NA NA NA -3.7 13.4 3.82
## 7 1990-01-10 NA NA NA NA NA 1.6 14.1 7.42
## 8 1990-01-11 NA NA NA NA NA -1.8 18.8 9.07
## 9 1990-01-12 NA NA NA NA NA -4.9 7.5 -0.0417
## 10 1990-01-13 NA NA NA NA NA -7.2 10.1 -0.229
## # ... with 10,944 more rows, and 13 more variables: RHMINcs21 <dbl>,
## # RHMAXcs21 <dbl>, RHcs21 <dbl>, WINDcs21 <dbl>, WDIRcs21 <dbl>,
## # TMINcs28 <dbl>, TMAXcs28 <dbl>, TAVGcs28 <dbl>, RHMINcs28 <dbl>,
## # RHMAXcs28 <dbl>, RHcs28 <dbl>, WINDcs28 <dbl>, WDIRcs28 <dbl>trimmedstartdaterhessys<- paste(year(ymd(trimmed_cwtclimate$date[1])), month(ymd(trimmed_cwtclimate$date[1])), day(ymd(trimmed_cwtclimate$date[1])),"1")
trimmedstartdaterhessys## [1] "1990 1 4 1"trimmed_cwtclimate <- trimmed_cwtclimate %>% select(-YEAR,-MONTH,-DAY)
# Replace NA values according to the specified rules
##NA APPROX is used, zoo package, interpolation
trimmed_cwtclimate <- trimmed_cwtclimate %>%
mutate(
RRG12 = ifelse(is.na(RRG12), 0, RRG12),
RRG20 = ifelse(is.na(RRG20), 0, RRG20),
TMINcs21 = ifelse(is.na(TMINcs21),TAVGcs21,TMINcs21),
TMINcs28 = ifelse(is.na(TMINcs28),TAVGcs28,TMINcs28),
TMINcs21 = na.approx(TMINcs21),
TMAXcs21 = na.approx(TMAXcs21),
TMINcs28 = na.approx(TMINcs28),
TMAXcs28 = na.approx(TMAXcs28),
TAVGcs21 = na.approx(TAVGcs21),
TAVGcs28 = na.approx(TAVGcs28),
TAVGcs21 = ifelse(TAVGcs21>TMINcs21 & TAVGcs21<TMAXcs21,TAVGcs21,(TMAXcs21+TMINcs21/2)),
TAVGcs28 = ifelse(TAVGcs28>TMINcs28 & TAVGcs28<TMAXcs28,TAVGcs28,(TMAXcs28+TMINcs28/2)),
RHMINcs21= na.approx(RHMINcs21),
RHMINcs28= na.approx(RHMINcs28),
RHMAXcs21= na.approx(RHMAXcs21),
RHMAXcs28= na.approx(RHMAXcs28),
RHcs21 = na.approx(RHcs21),
RHcs28 = na.approx(RHcs28),
WINDcs21 = na.approx(WINDcs21),
WINDcs28 = na.approx(WINDcs28),
WDIRcs21 = na.approx(WDIRcs21),
WDIRcs28 = na.approx(WDIRcs28)
)
# View the resulting dataframe
print(trimmed_cwtclimate)## # A tibble: 10,954 x 19
## date RRG12 RRG20 TMINcs21 TMAXcs21 TAVGcs21 RHMINcs21 RHMAXcs21
## <date> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1990-01-04 0.0429 0.0325 6.9 11.7 9.71 0.88 1
## 2 1990-01-05 0.00533 0.00381 5.5 14.3 8.73 0.7 1
## 3 1990-01-06 0.00686 0.00737 5.5 11.1 7.00 0.8 1
## 4 1990-01-07 0.0112 0.00991 5 6.8 5.9 0.91 1
## 5 1990-01-08 0.00584 0.00406 -0.2 13.2 6.46 0.68 1
## 6 1990-01-09 0 0 -3.7 13.4 3.82 0.41 1
## 7 1990-01-10 0 0 1.6 14.1 7.42 0.2 0.74
## 8 1990-01-11 0 0 -1.8 18.8 9.07 0.31 0.73
## 9 1990-01-12 0 0 -4.9 7.5 -0.0417 0.31 0.61
## 10 1990-01-13 0 0 -7.2 10.1 -0.229 0.22 0.6
## # ... with 10,944 more rows, and 11 more variables: RHcs21 <dbl>,
## # WINDcs21 <dbl>, WDIRcs21 <dbl>, TMINcs28 <dbl>, TMAXcs28 <dbl>,
## # TAVGcs28 <dbl>, RHMINcs28 <dbl>, RHMAXcs28 <dbl>, RHcs28 <dbl>,
## # WINDcs28 <dbl>, WDIRcs28 <dbl>rows_with_na <- trimmed_cwtclimate %>%
filter_all(any_vars(is.na(.)))
print(rows_with_na)## # A tibble: 0 x 19
## # ... with 19 variables: date <date>, RRG12 <dbl>, RRG20 <dbl>, TMINcs21 <dbl>,
## # TMAXcs21 <dbl>, TAVGcs21 <dbl>, RHMINcs21 <dbl>, RHMAXcs21 <dbl>,
## # RHcs21 <dbl>, WINDcs21 <dbl>, WDIRcs21 <dbl>, TMINcs28 <dbl>,
## # TMAXcs28 <dbl>, TAVGcs28 <dbl>, RHMINcs28 <dbl>, RHMAXcs28 <dbl>,
## # RHcs28 <dbl>, WINDcs28 <dbl>, WDIRcs28 <dbl>##make sure units are correct
cwtws32localrain<- data.frame(rbind(trimmedstartdaterhessys, as.data.frame(trimmed_cwtclimate$RRG12)))
cwtws32localtmin<- data.frame(rbind(trimmedstartdaterhessys, as.data.frame(trimmed_cwtclimate$TMINcs28)))
cwtws32localtmax<- data.frame(rbind(trimmedstartdaterhessys, as.data.frame(trimmed_cwtclimate$TMAXcs28)))
cwtws32localtavg <- data.frame(rbind(trimmedstartdaterhessys, as.data.frame(trimmed_cwtclimate$TAVGcs28)))
cwtws32localrhmin <- data.frame(rbind(trimmedstartdaterhessys, as.data.frame(trimmed_cwtclimate$RHMINcs28)))
cwtws32localrhmax <- data.frame(rbind(trimmedstartdaterhessys, as.data.frame(trimmed_cwtclimate$RHMAXcs28)))
cwtws32localrhavg <- data.frame(rbind(trimmedstartdaterhessys, as.data.frame(trimmed_cwtclimate$RHcs28)))
cwtws32localwind <- data.frame(rbind(trimmedstartdaterhessys, as.data.frame(trimmed_cwtclimate$WINDcs28)))
cwtws32localwdir <- data.frame(rbind(trimmedstartdaterhessys, as.data.frame(trimmed_cwtclimate$WDIRcs28)))
write.table(cwtws32localrain, file = "clim/cwtws32local.rain", quote = FALSE, row.names = FALSE, col.names = FALSE)
write.table(cwtws32localtmin, file = "clim/cwtws32local.tmin", quote = FALSE, row.names = FALSE, col.names = FALSE)
write.table(cwtws32localtmax, file = "clim/cwtws32local.tmax", quote = FALSE, row.names = FALSE, col.names = FALSE)
write.table(cwtws32localtavg, file = "clim/cwtws32local.tavg", quote = FALSE, row.names = FALSE, col.names = FALSE)
write.table(cwtws32localrhmin, file = "clim/cwtws32local.relative_humidity_min", quote = FALSE, row.names = FALSE, col.names = FALSE)
write.table(cwtws32localrhmax, file = "clim/cwtws32local.relative_humidity_max", quote = FALSE, row.names = FALSE, col.names = FALSE)
write.table(cwtws32localrhavg, file = "clim/cwtws32local.relativehumidity", quote = FALSE, row.names = FALSE, col.names = FALSE)
write.table(cwtws32localwind, file = "clim/cwtws32local.wind", quote = FALSE, row.names = FALSE, col.names = FALSE)
write.table(cwtws32localwdir, file = "clim/cwtws32local.wind_direction", quote = FALSE, row.names = FALSE, col.names = FALSE)Extend Local Datasets
setwd(system.file("extdata/", package = "RHESSysIOinR"))
precip <- read.csv('clim/cwtws32local.rain', header = FALSE)
tmin <- read.csv('clim/cwtws32local.tmin', header = FALSE)
tmax <- read.csv('clim/cwtws32local.tmax', header = FALSE)
head(precip)## V1
## 1 1990 1 4 1
## 2 0.042926
## 3 0.005334
## 4 0.006858
## 5 0.011176
## 6 0.005842nrow(precip)-1## [1] 10954substr(precip[1,1],1,nchar(precip[1,1])-2)## [1] "1990 1 4"dt <- as.Date(substr(precip[1,1],1,nchar(precip[1,1])-2), "%Y %m %d")
dt.new<- dt - as.difftime((nrow(precip)-1)*8, units="days")
#remove first value and repeats dataset 9 times
longcwtraindata<- do.call("rbind",replicate(9,as.list((precip[-1,]))))
longcwttmindata<- do.call("rbind",replicate(9,as.list((tmin[-1,]))))
longcwttmaxdata<- do.call("rbind",replicate(9,as.list((tmax[-1,]))))
longcwstartdate<- paste(year(ymd(dt.new)),month(ymd(dt.new)),day(ymd(dt.new)),"1")
longcwtrain<- rbind(longcwstartdate,longcwtraindata)
longcwttmin<- rbind(longcwstartdate,longcwttmindata)
longcwttmax<- rbind(longcwstartdate,longcwttmaxdata)
longcwtrain<- data.frame(longcwtrain)
longcwttmin<- data.frame(longcwttmin)
longcwtmax<- data.frame(longcwttmax)
write.table(longcwtrain, file = "clim/cwtws32extended.rain", quote = FALSE, row.names = FALSE, col.names = FALSE)
write.table(longcwttmin, file = "clim/cwtws32extended.tmin", quote = FALSE, row.names = FALSE, col.names = FALSE)
write.table(longcwttmax, file = "clim/cwtws32extended.tmax", quote = FALSE, row.names = FALSE, col.names = FALSE)
read_clim("clim/cwtws32extended.rain", dates_out=TRUE)## [1] "1750 01 30 01" "2019 12 31 24"read_clim("clim/cwtws32local.rain", dates_out=TRUE)## [1] "1990 01 04 01" "2019 12 31 24"Create DAYMET Base File
setwd(system.file("extdata/", package = "RHESSysIOinR"))
##Quick way to create base file necessary for climate processing
{daymetbasefiletext<- "101 base_station_id
275535 x_coordinate
3879165 y_coordinate
1064 z_coordinate
3.5 effective_lai
5.0 screen_height
clim/cwt_annual annual_climate_prefix
0 number_non_critical_annual_sequences
clim/cwt_monthly monthly_climate_prefix
0 number_non_critical_monthly_sequences
clim/cwtws32daymet daily_climate_prefix
0 number_non_critical_daily_sequences
clim/cwt_hourly hourly_climate_prefix
0 number_non_critical_hourly_sequences
clim/cwt trigger_of_dated_input
0 num_non_critical_sequences
"}
write(daymetbasefiletext,file = "clim/cwtws32daymet.base")Download Daymet Data for Watershed of Interest using daymetr package Currently only creates precip, tmin and tmax files.
Extend DAYMET datasets to allow for longer spinup period Read Daymet Datasets and Replicate 5x times to create a longer period of record for spin up Create base File manually for now
setwd(system.file("extdata/", package = "RHESSysIOinR"))
getwd()## [1] "C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata"precip <- read.csv('clim/cwtws32daymet.rain', header = FALSE)
tmin <- read.csv('clim/cwtws32daymet.tmin', header = FALSE)
tmax <- read.csv('clim/cwtws32daymet.tmax', header = FALSE)
#tavg <- read.csv('clim/cwtws32daymet.tavg', header = FALSE)
#rh <- read.csv('clim/cwtws32daymet.relative_humidity', header = FALSE)
#wind <- read.csv('clim/cwtws32daymet.wind', header = FALSE)
head(precip)## V1
## 1 1980 1 1 1
## 2 0.00427
## 3 0
## 4 0.00866
## 5 0.01468
## 6 0nrow(precip)-1## [1] 15341substr(precip[1,1],1,nchar(precip[1,1])-2)## [1] "1980 1 1"dt <- as.Date(substr(precip[1,1],1,nchar(precip[1,1])-2), "%Y %m %d")
dt.new<- dt - as.difftime((nrow(precip)-1)*5, units="days")
#remove first value and repeats dataset 6 times
longcwtraindata<- do.call("rbind",replicate(6,as.list((precip[-1,]))))
longcwttmindata<- do.call("rbind",replicate(6,as.list((tmin[-1,]))))
longcwttmaxdata<- do.call("rbind",replicate(6,as.list((tmax[-1,]))))
#longcwttavgdata<- do.call("rbind",replicate(6,as.list((tavg[-1,]))))
#longcwtrhdata<- do.call("rbind",replicate(6,as.list((rh[-1,]))))
#longcwtwinddata<- do.call("rbind",replicate(6,as.list((wind[-1,]))))
longcwstartdate<- paste(year(ymd(dt.new)),month(ymd(dt.new)),day(ymd(dt.new)),"1")
longcwtrain<- rbind(longcwstartdate,longcwtraindata)
longcwttmin<- rbind(longcwstartdate,longcwttmindata)
longcwttmax<- rbind(longcwstartdate,longcwttmaxdata)
#longcwttavg<- rbind(longcwstartdate,longcwttavgdata)
#longcwtrh<- rbind(longcwstartdate,longcwtrhdata)
#longcwtwind<- rbind(longcwstartdate,longcwtwinddata)
longcwtrain<- data.frame(longcwtrain)
longcwttmin<- data.frame(longcwttmin)
longcwtmax<- data.frame(longcwttmax)
#longcwttavg<- data.frame(longcwttavg)
#longcwtrh<- data.frame(longcwtrh)
#longcwtwind<- data.frame(longcwtwind)
write.table(longcwtrain, file = "clim/cwtws32daymetextended.rain", quote = FALSE, row.names = FALSE, col.names = FALSE)
write.table(longcwttmin, file = "clim/cwtws32daymetextended.tmin", quote = FALSE, row.names = FALSE, col.names = FALSE)
write.table(longcwttmax, file = "clim/cwtws32daymetextended.tmax", quote = FALSE, row.names = FALSE, col.names = FALSE)
#write.table(longcwttavg, file = "clim/cwtws32extended.tavg", quote = FALSE, row.names = FALSE, col.names = FALSE)
#write.table(longcwtrh, file = "clim/cwtws32extended.relative_humidity", quote = FALSE, row.names = FALSE, col.names = FALSE)
#write.table(longcwtwind, file = "clim/cwtws32extended.wind", quote = FALSE, row.names = FALSE, col.names = FALSE)
read_clim("clim/cwtws32extended.rain", dates_out=TRUE)## [1] "1750 01 30 01" "2019 12 31 24"read_clim("clim/cwtws32daymet.rain", dates_out=TRUE)## [1] "1980 01 01 01" "2021 12 31 24"Download and Prepare DEM Code to download DEM originally prepared by Dylan Beaudette
setwd(system.file("extdata/", package = "RHESSysIOinR"))
# use watershed polygons for BBOX to request elevation data
x <- read_sf(watershedshapefile)
# buffer to ensure that there are no truncated watersheds
x.buff <- st_as_sf(st_buffer(st_as_sfc(st_bbox(x)), dist = 1000))
# convert to GCS WGS84 -> no transformation / warp will be applied to DEM
x.gcs <- st_transform(x, 4326)
x.buff.gcs <- st_transform(x.buff, 4326)
# requires sf collection (geometry + attributes)
# get DEM in GCS WGS84 and use z = 14 for best available data
elevationraster <- get_elev_raster(locations = x.buff.gcs, z = 14, clip = 'bbox')
res(elevationraster)## [1] 3.921794e-05 3.921794e-05# check: OK
{plot(elevationraster, main = "DEM, Watershed Outline and Buffer" )
plot(st_geometry(x.buff.gcs), add = TRUE, col = NA, lwd = 2)
plot(st_geometry(x.gcs), add = TRUE, col = NA)}
{plot(elevationraster, main = "DEM and selected Watershed")
plot(st_geometry(x.gcs[which(x$WS=='32'),]), add = TRUE, col = NA)}
cropped_dem<- crop(elevationraster,st_bbox(x.gcs))
plot(cropped_dem, main = "DEM cropped to Buffer")
plot(st_geometry(x.gcs[which(x$WS=='32'),]), add = TRUE, col = NA)
newproj<- "+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs"
reproje<- projectRaster(from = cropped_dem, crs = newproj, method = 'bilinear')
plot(reproje, main = "Reprojected cropped DEM")
dummydata<- reproje
res(dummydata) <- 10
resampleddem <- resample(reproje, dummydata, method = 'bilinear')
plot(resampleddem, main = "DEM reprojected to 10m Resolution")
{plot(resampleddem, main = "DEM reprojected to 10m Resolution")
plot(st_geometry(x.buff.gcs), add = TRUE, col = NA, lwd = 2)
plot(st_geometry(x.gcs), add = TRUE, col = NA)}
Start Spatial Processing in GRASS8 Code to prepare spatial input files for RHESSYs takes cues from Will Burkes “spatial_input_gen.R” script but is prepared for GRASS8 instead of GRASS7
setwd(system.file("extdata/", package = "RHESSysIOinR"))
#gisbase_grass <- ifelse(.Platform$OS.type == "windows", link2GI::paramGRASSw()$gisbase_GRASS[1], link2GI::paramGRASSx()$gisbase_GRASS[1])
## set easting and northing as coordinates of lowest point within boundary basin, set threshold for basin delineation
threshold = 1000
easting = NULL
northing = NULL
dir.create("grassdata")
initGRASS(gisBase = "C:/Program Files/GRASS GIS 8.2",
#gisDbase = "C:/Users/Carlos/Documents/grassdata",
gisDbase = "grassdata",
location = "Coweeta1",
mapset = "PERMANENT",
override = TRUE)## gisdbase grassdata
## location Coweeta1
## mapset PERMANENT
## rows 611
## columns 537
## north 3884060
## south 3877950
## west 273846.7
## east 279216.7
## nsres 10
## ewres 10
## projection +proj=utm +no_defs +zone=17 +a=6378137 +rf=298.257223563
## +towgs84=0.000,0.000,0.000 +type=crs +to_meter=1execGRASS("g.proj", proj4="+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs", flag = "c" )
execGRASS("g.gisenv", flag="n")
execGRASS("g.region", flag=c("p"))
##Convert DEM to Spatrast and then import DEM to GRASS
spatrast_dem10 <- rast(resampleddem)
write_RAST(spatrast_dem10, "grass_dem10", overwrite = TRUE)
## The following code is heavily influenced by Will Burkes "spatial_input_gen.R" script
execGRASS("g.region", raster = "grass_dem10", flag=c("p"))
## Import shapefiles (Streams, Weirs, Watersheds, Roads)
## only roads and Watersheds necessary to Run
streams<- vect(streamshapefile)
streams<- project(streams,"+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs")
write_VECT(streams, "streams")
weirs<- vect(weirshapefile)
weirs<- project(weirs,"+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs")
write_VECT(weirs, "weirs")
watersheds<- vect(watershedshapefile)
wgs84watersheds<- project(watersheds,"+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs")
write_VECT(wgs84watersheds, "watersheds")
roads <- vect(roadshapefile)
roads<- project(roads,"+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs")
write_VECT(roads,"roads")
## Check all vectors that have been exported to GRASS
execGRASS("g.list", parameters = list(type = "vector"))
## Generate a depressionless DEM and d-8 flow direction maps from DEM
execGRASS("r.fill.dir", input="grass_dem10", output="dem10f", direction="dir10", flags = c("overwrite") )
## Generate flow accumulation, drainage, and horizon maps
execGRASS("r.watershed", elevation = "dem10f", accumulation = "acc10", drainage = "drain10")
## Generate slope, aspect, wetness index, horizon, and road maps
execGRASS("r.slope.aspect", elevation = "grass_dem10", slope = "slope10", aspect = "aspect10")
execGRASS("r.topidx", input = "grass_dem10", output = "topidx10")
execGRASS("r.mapcalc", expression = "topidx10_100 = topidx10 * 100", flags = "overwrite")
execGRASS("r.horizon", flags = "d", elevation = "grass_dem10", direction = 0, output = "east")
execGRASS("r.horizon", flags = "d", elevation = "grass_dem10", direction = 180, output = "west")
execGRASS("r.mapcalc", expression = "e_horizon_1000 = sin(east_000) * 1000", flags = "overwrite")
execGRASS("r.mapcalc", expression = "w_horizon_1000 = sin(west_180) * 1000", flags = "overwrite")
## Generate a hillslope and stream map
subbasin_name = paste0("subbasin_th",as.character(threshold))
stream_name = paste0("stream_th",as.character(threshold))
hillslope_name = paste0("hillslope_th",as.character(threshold))
execGRASS("r.watershed", elevation = "dem10f", threshold = threshold, basin = subbasin_name, stream = stream_name, half_basin = hillslope_name)
## Delineate watershed boundary
weirs[12,]## class : SpatVector
## geometry : points
## dimensions : 1, 1 (geometries, attributes)
## extent : 276120.4, 276120.4, 3881339, 3881339 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=17 +datum=WGS84 +units=m +no_defs
## names : id
## type : <int>
## values : 12#easting = crds(weirs[12,])[1]
#northing = crds(weirs[12,])[2]
# should snap to stream
easting = 276111
northing = 3881335
execGRASS("g.region", raster = "grass_dem10", flag=c("p"))
execGRASS("r.water.outlet", input ="drain10", output="basin_ws32",coordinates = c(easting, northing), flags = "overwrite")
## Mask Data
# execGRASS("r.mask", input = "basin_ws32", maskcats = 1)
execGRASS("r.info",map = "basin_ws32")
execGRASS("r.info",map = "grass_dem10")
## Generate Patch Maps
execGRASS("r.mapcalc", expression = "patch = row() * 537 + col()", flags = "overwrite")
## Create Road Maps
execGRASS("v.to.rast", input="roads" ,output="roads", use="cat", label_column="ROADNAME", flags = "overwrite")
execGRASS("r.mapcalc", expression = "roads = roads > 0", flags = "overwrite")
execGRASS("r.null", map = "roads", null= 0)
## Check all rasters that have been exported to or created in GRASS
execGRASS("g.list", parameters = list(type = "raster"))
## Export Prepared Maps from GRASS to R
filleddem <- read_RAST("dem10f")
##Export maps back into R or export directly from GRASS to disk
grass_acc10 <- read_RAST("acc10")
grass_aspect10 <- read_RAST("aspect10")
grass_basin_ws32 <- read_RAST("basin_ws32")
grass_dem10 <- read_RAST("grass_dem10")
grass_dem10f <- read_RAST("dem10f")
grass_dir10 <- read_RAST("dir10")
grass_drain10 <- read_RAST("drain10")
grass_e_horizon_1000 <- read_RAST("e_horizon_1000")
grass_hillslope_th1000 <- read_RAST("hillslope_th1000")
#grass_patch <- read_RAST("patch")
grass_slope10 <- read_RAST("slope10")
grass_stream_th1000 <- read_RAST("stream_th1000")
grass_topidx10 <- read_RAST("topidx10")
grass_topidx10_100 <- read_RAST("topidx10_100")
grass_w_horizon_1000 <- read_RAST("w_horizon_1000")
grass_east_000 <- read_RAST("east_000")
grass_west_180 <- read_RAST("west_180")
grass_subbasin_th1000 <- read_RAST("subbasin_th1000")
dir.create("grassexport")
dir.create("grassexport/cwt")
execGRASS("r.out.gdal", input = "acc10", output = "grassexport/cwt/acc10.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "aspect10", output = "grassexport/cwt/aspect10.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "basin_ws32", output = "grassexport/cwt/basin_ws32.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "grass_dem10", output = "grassexport/cwt/grass_dem10.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "dem10f", output = "grassexport/cwt/dem10f.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "dir10", output = "grassexport/cwt/dir10.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "drain10", output = "grassexport/cwt/drain10.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "e_horizon_1000", output = "grassexport/cwt/e_horizon_1000.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "hillslope_th1000", output = "grassexport/cwt/hillslope_th1000.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "patch", output = "grassexport/cwt/patch.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "slope10", output = "grassexport/cwt/slope10.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "stream_th1000", output = "grassexport/cwt/stream_th1000.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "topidx10", output = "grassexport/cwt/topidx10.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "topidx10_100", output = "grassexport/cwt/topidx10_100.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "w_horizon_1000", output = "grassexport/cwt/w_horizon_1000.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "east_000", output = "grassexport/cwt/east_000.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "west_180", output = "grassexport/cwt/west_180.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "subbasin_th1000", output = "grassexport/cwt/subbasin_th1000.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "roads", output = "grassexport/cwt/roads.tif", format = "GTiff", flags = "overwrite")Crop Maps that were prepared in GRASS. Maps are cropped R because MASK command is not functioning properly for GRASS via R. Original input maps are been cropped before processing worldfile and flowtable in R, without cropping the RHESSYs worldfile and flowtable preparation in the R environment will not function properly,
setwd(system.file("extdata/", package = "RHESSysIOinR"))
ws32<- raster("grassexport/cwt/basin_ws32.tif")
acc10<- raster("grassexport/cwt/acc10.tif")
patches<- raster("grassexport/cwt/patch.tif")
aspect10<- raster("grassexport/cwt/aspect10.tif")
dem10<- raster("grassexport/cwt/grass_dem10.tif")
dem10f<- raster("grassexport/cwt/dem10f.tif")
dir10<- raster("grassexport/cwt/dir10.tif")
drain10<- raster("grassexport/cwt/drain10.tif")
hillslope<- raster("grassexport/cwt/hillslope_th1000.tif")
roads<- brick("grassexport/cwt/roads.tif")
slope10<- raster("grassexport/cwt/slope10.tif")
streams<- raster("grassexport/cwt/stream_th1000.tif")
topidx10_100<- raster("grassexport/cwt/topidx10_100.tif")
e_horizon_1000 <- raster("grassexport/cwt/e_horizon_1000.tif")
w_horizon_1000 <- raster("grassexport/cwt/w_horizon_1000.tif")
subbasins<- raster("grassexport/cwt/subbasin_th1000.tif")
ws32polygon<- rasterToPolygons(ws32, fun = function(x){x==1}, dissolve = TRUE, digits = 12)
maskedws32 <- crop(mask(ws32,ws32polygon),ws32polygon)
maskedacc10 <- crop(mask(acc10,ws32polygon),ws32polygon)
maskedpatches <- crop(mask(patches,ws32polygon),ws32polygon)
maskedaspect10 <- crop(mask(aspect10,ws32polygon),ws32polygon)
maskeddem10 <- crop(mask(dem10,ws32polygon),ws32polygon)
maskeddem10f <- crop(mask(dem10f,ws32polygon),ws32polygon)
maskeddir10 <- crop(mask(dir10,ws32polygon),ws32polygon)
maskeddrain10 <- crop(mask(drain10,ws32polygon),ws32polygon)
maskedhillslope <- crop(mask(hillslope,ws32polygon),ws32polygon)
maskedroads <- crop(mask(roads,ws32polygon),ws32polygon)
maskedslope10 <- crop(mask(slope10,ws32polygon),ws32polygon)
maskedstreams <- crop(mask(streams,ws32polygon),ws32polygon)
maskedtopidx10_100 <- crop(mask(topidx10_100,ws32polygon),ws32polygon)
maskede_horizon_1000 <- crop(mask(e_horizon_1000,ws32polygon),ws32polygon)
maskedw_horizon_1000 <- crop(mask(w_horizon_1000,ws32polygon),ws32polygon)
maskedsubbasins<- crop(mask(subbasins,ws32polygon),ws32polygon)
## Save as tif
dir.create("grassexport/cwt_ws32")
writeRaster(maskedws32,"grassexport/cwt_ws32/basin_ws32",format = "GTiff", overwrite = TRUE)
writeRaster(maskedacc10,"grassexport/cwt_ws32/acc10",format = "GTiff", overwrite = TRUE)
writeRaster(maskedpatches,"grassexport/cwt_ws32/patch",format = "GTiff", overwrite = TRUE)
writeRaster(maskedaspect10,"grassexport/cwt_ws32/aspect10",format = "GTiff", overwrite = TRUE)
writeRaster(maskeddem10,"grassexport/cwt_ws32/dem10",format = "GTiff", overwrite = TRUE)
writeRaster(maskeddem10f,"grassexport/cwt_ws32/dem10f",format = "GTiff", overwrite = TRUE)
writeRaster(maskeddir10,"grassexport/cwt_ws32/dir10",format = "GTiff", overwrite = TRUE)
writeRaster(maskeddrain10,"grassexport/cwt_ws32/drain10",format = "GTiff", overwrite = TRUE)
writeRaster(maskedhillslope,"grassexport/cwt_ws32/hillslope",format = "GTiff", overwrite = TRUE)
writeRaster(maskedroads,"grassexport/cwt_ws32/roads",format = "GTiff", overwrite = TRUE)
writeRaster(maskedslope10,"grassexport/cwt_ws32/slope10",format = "GTiff", overwrite = TRUE)
writeRaster(maskedstreams,"grassexport/cwt_ws32/streams",format = "GTiff", overwrite = TRUE)
writeRaster(maskedtopidx10_100,"grassexport/cwt_ws32/topidx10_100",format = "GTiff", overwrite = TRUE)
writeRaster(maskede_horizon_1000,"grassexport/cwt_ws32/e_horizon_1000",format = "GTiff", overwrite = TRUE)
writeRaster(maskedw_horizon_1000,"grassexport/cwt_ws32/w_horizon_1000",format = "GTiff", overwrite = TRUE)
writeRaster(maskedsubbasins,"grassexport/cwt_ws32/subbasins",format = "GTiff", overwrite = TRUE)
plot(maskedws32, main ="ws32")
plot(maskedpatches, main = "patches")
plot(maskedacc10, main = "acc10")
plot(maskedaspect10, main = "aspect10")
plot(maskedsubbasins, main = "subbasins")
plot(maskeddem10, main = "dem10")
plot(maskeddem10f, main = "dem10f")
plot(maskeddir10, main = "dir10")
plot(maskeddrain10, main = "drain10")
plot(maskedhillslope, main = "hillslope")
plot(maskedroads, main = "roads")
plot(maskedslope10, main = "slope10")
plot(maskedstreams, main = "streams")
plot(maskedtopidx10_100, main = "topidx10_100")
plot(maskede_horizon_1000, main = "e_horizon_1000")
plot(maskedw_horizon_1000, main = "w_horizon_1000")
# This must be properly intergrated into code # RSS files from github or folder? Prepare SSURGO Soil Maps SSURGO-via-SDA Code originally prepared by Dylan Beaudette Download SSURGO data and make sure CRS/Extent/Crop are Correct for input into RHESSys
setwd(system.file("extdata/", package = "RHESSysIOinR"))
# load BBOX
x <- read_sf(watershedshapefile)
# get gSSURGO grid here
mu <- mukey.wcs(aoi = x, db = 'gssurgo')
# unique map unit keys
ll <- levels(mu)[[1]]
# note SSA boundary
levelplot(mu, att = 'ID', margin = FALSE, colorkey = FALSE, col.regions = viridis, main = "SSURGO Map Units")
## Note: some of these map units are dominated by non-soil components
# Dominant Component (Numeric) -> NODATA
# Weighted Average -> use any available data, but wt. mean are diluted (see source data)
# get thematic data from SDA
# dominant component
# depth-weighted average
# sand, silt, clay, AWC (RV)
p <- get_SDA_property(property = c("sandtotal_r","silttotal_r","claytotal_r", "awc_r", "ksat_r"),
method = "Weighted Average",
mukeys = ll$ID,
top_depth = 0,
bottom_depth = 200)
head(p)## areasymbol musym
## 1 NC113 BuD
## 2 NC113 BuF
## 3 NC113 CaF
## 4 NC113 CdD
## 5 NC113 CdE
## 6 NC113 CdF
## muname
## 1 Burton-Craggey-Rock outcrop complex, windswept, 15 to 30 percent slopes, stony
## 2 Burton-Craggey-Rock outcrop complex, windswept, 30 to 95 percent slopes, stony
## 3 Cashiers gravelly fine sandy loam, 50 to 95 percent slopes
## 4 Chandler gravelly fine sandy loam, 15 to 30 percent slopes
## 5 Chandler gravelly fine sandy loam, 30 to 50 percent slopes
## 6 Chandler gravelly fine sandy loam, 50 to 95 percent slopes
## mukey sandtotal_r silttotal_r claytotal_r awc_r ksat_r
## 1 545800 66.36 19.50 14.14 0.05 9.78
## 2 545801 66.41 19.28 14.32 0.05 9.58
## 3 545803 66.82 21.68 11.50 0.14 28.00
## 4 545805 46.74 41.76 11.50 0.13 28.00
## 5 545806 46.74 41.76 11.50 0.13 28.00
## 6 545807 46.74 41.76 11.50 0.13 28.00# re-create raster attribute table with aggregate soil properties
rat <- merge(ll, p, by.x = 'ID', by.y = 'mukey', sort = FALSE, all.x = TRUE)
# re-pack RAT
levels(mu) <- rat
# convert raster + RAT --> stack of values
s <- deratify(mu, att = c("sandtotal_r", "silttotal_r", "claytotal_r", "awc_r", "ksat_r"))
# graphical check
levelplot(
s[[1:3]],
main = 'Sand, Silt, Clay (RV) 0-200cm\nWeighted Average',
margin = FALSE,
scales = list(draw = FALSE),
col.regions = viridis
)
levelplot(
s[[4]],
main = 'AWC (RV) 0-200cm\nWeighted Average',
margin = FALSE,
scales = list(draw = FALSE),
col.regions = viridis
)
levelplot(
s[[5]],
main = 'ksat (RV) 0-200cm\nWeighted Average',
margin = FALSE,
scales = list(draw = FALSE),
col.regions = viridis
)
# convert to a representative soil texture class
txt.lut <- read.csv('http://soilmap2-1.lawr.ucdavis.edu/800m_grids/RAT/texture_2550.csv')
# make a copy
texture_060 <- s[[1]]
# note: soil textures that aren't present are dropped from factor levels
texture_060[] <- ssc_to_texcl(sand = s$sandtotal_r[], clay = s$claytotal_r[])
# extract RAT
rat <- levels(texture_060)[[1]]
# add colors
rat <- merge(rat, txt.lut, by.x = 'VALUE', by.y = 'class', all.x = TRUE, sort = FALSE)
# fix column order
rat <- rat[, c('ID', 'VALUE', 'hex', 'names')]
# re-pack
levels(texture_060) <- rat
# check: ok
cols <- levels(texture_060)[[1]]$hex
levelplot(texture_060, att = 'names', margin = FALSE, col.regions = cols, main = "Soil Textures") 
dir.create("rasters/")
dir.create("rasters/cwt")
dir.create("rasters/cwt/ws32")
writeRaster(texture_060, file = 'rasters/cwt/ws32/soiltexture060.tif', overwrite = TRUE, options=c("COMPRESS=LZW"))
texture60raster<- raster("rasters/cwt/ws32/soiltexture060.tif")
### this has been modified to use new SSURGO Inputs with deep-shallow depth classes
texture60raster<- raster("rasters/cwt/ws32/ssurgo-soiltype-class.tif")
## reproject to raster mask instead
ws32<- raster("grassexport/cwt_ws32/basin_ws32.tif")
newproj<- "+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs"
reproj60 <- projectRaster(from = texture60raster, crs = proj4string(ws32), method = 'ngb')
plot(texture60raster, zlim = c(1,5), main = "Texture Raster with modified Depth Classes")
{plot(reproj60, zlim = c(1,5), main = "Reprojected Texture Raster with modified Depth Classes and Watersheds")
plot(st_geometry(x), add = TRUE, col = NA)}
## prepare 10m resolution raster and resample prepared soil data to 10m
resample10<- raster(resolution=c(10,10), crs = proj4string(reproj60), ext=extent(reproj60))
reproj6010 <- resample(reproj60,resample10, method = 'ngb')
{plot(reproj6010, zlim = c(1,5), main = "Reprojected and Resampled Texture Raster with modified Depth Classes and Watersheds")
plot(st_geometry(x), add = TRUE, col = NA)}
## Crop prepared raster to extent of watershed of interest
## croppedtexture60 <- crop(reproj6010,extent(x[which(x$WS=='32'),]))
#for now cropping is completed based on delineated watershed extent from raster instead of watershed extent from watershed polygons
croppedtexture60 <- crop(reproj6010,extent(ws32))
# croppedtexture60 <- projectRaster(from = croppedtexture60, crs = newproj, method = 'ngb')
## Mask prepared raster using watershed polygon
#masked60 <- mask(croppedtexture60,x[which(x$WS=='32'),])
## for now masking is completed based on delineated watershed extent from raster instead of watershed extent from watershed polygons
ws32polygon<- rasterToPolygons(ws32, fun = function(x){x==1}, dissolve = TRUE, digits = 12)
extent(croppedtexture60)<- extent(ws32)
{plot(reproj6010, zlim = c(1,5), main = "Reprojected and Resampled Texture Raster with modified Depth Classes and Selected Watershed")
plot(ws32polygon, add = TRUE, col = NA)}
masked60 <- mask(croppedtexture60,ws32polygon)
plot(masked60, zlim = c(1,5), main = "Reprojected and Resampled Texture Raster with modified Depth Classes and Masked Selected Watershed")
plot(x[which(x$WS=='32'),])
plot(st_geometry(x[which(x$WS=='32'),]))
{plot(reproj60, zlim = c(1,5))
plot(st_geometry(x[which(x$WS=='32'),]), add = TRUE, col = NA)}
{plot(as.factor(croppedtexture60), zlim = c(1,5))
plot(st_geometry(x[which(x$WS=='32'),]), add = TRUE, col = NA)}
{plot(as.factor(masked60), zlim = c(1,5))
plot(st_geometry(x[which(x$WS=='32'),]), add = TRUE, col = NA)}
{plot(masked60, zlim = c(1,5), main = "Polygon Delineation")
plot(st_geometry(x[which(x$WS=='32'),]), add = TRUE, col = NA)}
{plot(masked60, zlim = c(1,5), main = "Watershed Delineation")
plot(ws32polygon, add = TRUE, col = NA)}
Prepare textural classes from sand and clay totals
mapunittable<- levels(mu)[[1]]
paged_table(mapunittable)mapunittable$textclass<- ssc_to_texcl(sand = mapunittable$sandtotal_r, clay = mapunittable$claytotal_r)
paged_table(mapunittable)Code Prepared By Dylan Beaudette and modified
setwd(system.file("extdata/", package = "RHESSysIOinR"))
ws32<- raster("grassexport/cwt_ws32/basin_ws32.tif")
#("grassexport/cwt_ws32/ws32.tif")
# RSS: cell values are map unit keys
r <- raster('Raster soil survey/Coweeta_Final_raster.tif')
## for now, convert to UTM z17
r <- projectRaster(r, ws32, method = 'ngb')
# convert to raster + RAT
r <- ratify(r)
# RAT: raster attribute table
rat <- read.dbf('Raster soil survey/Coweeta_Final_raster.tif.vat.dbf', as.is = TRUE)
names(rat) <- tolower(names(rat))
# musym is the national mu symbol
rss.mu <- st_read('Raster soil survey/rss_nc/rss_nc.gdb', layer = 'mapunit')## Reading layer `mapunit' from data source
## `C:\Users\Carlos\Documents\R\win-library\4.0\RHESSysIOinR\extdata\Raster Soil Survey\rss_nc\rss_nc.gdb'
## using driver `OpenFileGDB'rss.co <- st_read('Raster soil survey/rss_nc/rss_nc.gdb', layer = 'component')## Reading layer `component' from data source
## `C:\Users\Carlos\Documents\R\win-library\4.0\RHESSysIOinR\extdata\Raster Soil Survey\rss_nc\rss_nc.gdb'
## using driver `OpenFileGDB'rss.hz <- st_read('Raster soil survey/rss_nc/rss_nc.gdb', layer = 'chorizon')## Reading layer `chorizon' from data source
## `C:\Users\Carlos\Documents\R\win-library\4.0\RHESSysIOinR\extdata\Raster Soil Survey\rss_nc\rss_nc.gdb'
## using driver `OpenFileGDB'## check for missing symbols
# none missing from RAT
setdiff(rat$MUSYM, rss.mu$musym)## NULL# map unit table contains some extra
setdiff(rss.mu$musym, rat$MUSYM)## [1] "2y9mf" "2y9mh" "2y9mj" "2y9mk" "2y9ml" "2y9mm" "2y9mn" "2y9mp" "2y9mq"
## [10] "2y9mr" "2y9ms" "2y9mt" "2y9mv" "2y9mw" "2y9mx" "2y9mg" "2y9pb" "2y9pc"
## [19] "2y9pd" "2y9p9" "2y9nz" "2y9nl" "2y9nm" "2y9nw" "2y9nx" "2zyxq" "2zyxw"
## [28] "2zyxy" "2zyy5" "2zyxz" "2zyy0" "2zyxx" "2zyxt" "2zyy6" "2y9np" "2y9nt"
## [37] "2y9nv" "2zyxr" "2zyy1" "2zyy2" "2zyy3" "2y9nh" "2y9nj" "2y9nk" "2zyxv"
## [46] "2zyyb" "2zyyc" "2zyyd" "2zyyf" "2y9nn" "2y9p1" "2y9p2" "2y9nq" "2y9nr"
## [55] "2zyy4"## subset child tables
rss.mu <- rss.mu[rss.mu$musym %in% rat$musym, ]
rss.co <- rss.co[rss.co$mukey %in% rss.mu$mukey, ]
rss.hz <- rss.hz[rss.hz$cokey %in% rss.co$cokey, ]
.dominantCondition <- function(i, v) {
i <- i[i$compkind != 'Miscellaneous area', ]
if(nrow(i) < 1) {
return(NULL)
}
fm <- as.formula(sprintf("comppct_r ~ %s", v))
a <- aggregate(fm, data = i, FUN = sum, na.rm = TRUE)
idx <- order(a[['comppct_r']], decreasing = TRUE)[1]
res <- data.frame(
mukey = i$mukey[1],
v = a[[v]][idx]
)
names(res) <- c('mukey', v)
return(res)
}
.dominantValue <- function(i, v) {
i <- i[i$compkind != 'Miscellaneous area', ]
if(nrow(i) < 1) {
return(NULL)
}
idx <- order(i[['comppct_r']], decreasing = TRUE)[1]
res <- data.frame(
mukey = i$mukey[1],
v = i[[v]][idx]
)
names(res) <- c('mukey', v)
return(res)
}
getDominantCondition <- function(x, v) {
res <- lapply(
split(rss.co, rss.co$mukey),
.dominantCondition, v = v
)
res <- do.call('rbind', res)
return(res)
}
getDominantValue <- function(x, v) {
res <- lapply(
split(rss.co.aws050, rss.co.aws050$mukey),
.dominantValue, v = v
)
res <- do.call('rbind', res)
return(res)
}
co.taxpartsize <- getDominantCondition(rss.co, v = 'taxpartsize')
co.spc <- rss.hz
depths(co.spc) <- cokey ~ hzdept_r + hzdepb_r
hzdesgnname(co.spc) <- 'hzname'
# at dZ = 1cm, use awc_r directly
co.spc$aws <- co.spc$awc_r * 1
a <- slab(co.spc, cokey ~ aws, slab.structure = c(0, 50), slab.fun = sum)
head(a)## variable cokey value top bottom contributing_fraction
## 1 aws 3244721:2807027 2.40 0 50 1
## 2 aws 3244721:2807439 4.02 0 50 1
## 3 aws 3244721:2807442 2.50 0 50 1
## 4 aws 3244721:2807445 3.73 0 50 1
## 5 aws 3244722:2760595 2.40 0 50 1
## 6 aws 3244722:2760596 4.02 0 50 1rss.co.aws050 <- merge(rss.co, a, by = 'cokey', sort = FALSE)[, c('mukey', 'cokey', 'compkind', 'comppct_r', 'value')]
# no NA allowed in wt. mean / etc.
rss.co.aws050 <- rss.co.aws050[!is.na(rss.co.aws050$value), ]
co.aws050 <- getDominantValue(rss.co.aws050, v = 'value')
names(co.aws050) <- c('mukey', 'aws050')
mu.subset <- rss.mu[, c('mukey', 'musym', 'mukind')]
agg <- merge(mu.subset, co.taxpartsize, by = 'mukey', sort = FALSE)
agg <- merge(agg, co.aws050, by = 'mukey', sort = FALSE)
rat <- merge(rat, agg, by = 'musym', sort = FALSE)
# attempt to split out component / series names
rat$co.names <- stri_split_fixed(str = rat$muname, pattern = ',', n = 2, simplify = TRUE)[, 1]
# this only works because all component names are single-word names
rat$co.names <- stri_split_fixed(str = rat$co.names, pattern = ' ', n = 2, simplify = TRUE)[, 1]
# RAT management
ll.original <- levels(r)[[1]]
# merge and pack updated RAT
ll <- merge(ll.original, rat, by.x = 'ID', by.y = 'value', all.x = TRUE, sort = FALSE)
levels(r) <- ll
# convert select attributes to new raster objects via RAT + codes
r.taxpartsize <- deratify(r, att = 'taxpartsize')
plot(r.taxpartsize)
###quick and dirty incorporate new DSM map with two depth classes
r.taxpartsizetest<- raster("C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata/Raster Soil Survey/rss-soiltype-class.tif")
r.taxpartsize <- projectRaster(r.taxpartsizetest, ws32, method = 'ngb')
ws32polygon<- rasterToPolygons(ws32, fun = function(x){x==1}, dissolve = TRUE, digits = 12)
maskeddsm <- mask(r.taxpartsize,ws32polygon)
plot(maskeddsm, main = "Masked DSM")
Reclassify Soil Maps
Default Texture Values in RHESSys 1 - Clay 2 - Silty Clay 3 - Silty Clay Loam 4 - Sandy Clay 5 - Sandy Clay Loam 6 - Clay Loam 7 - Silt 8 - Silty Loam 9 - Loam 10 - Sand 11 - Loamy Sand 12 - Sandy Loam 13 - Rock 14 - Water
Shallow Files Created and added to defs folder 19 - Shallow Sandy Clay Loam 23 - Shallow Loam 26 - Shallow Sandy Loam
Default Texture Values used by SDA 1 - loamy sand 2 - sandy loam 3 - loam 4 - sandy clay loam 5 - clay loam
Default Texture used by DSM and quick and dirty conversions 1 - loamy > loam 2 - fine-loamy > silty clay loam 3 - coarse-loamy > sandy loam
Current Textural Levels for Coweeta 6 Levels 1 - deep.sl 2- shallow.sl 3- deep.l 4- shallow.l 5- deep.scl 6- shallow.scl
setwd(system.file("extdata/", package = "RHESSysIOinR"))
# original before depth class substitution
#reclasssoil<- c(1,11,2,12,3,9,4,5,5,6,6,0,7,0,8,0,9,0,10,0,11,0,12,0,13,0,14,0)
##this has been modified for quick implementation of shallow-deep classes
reclasssoil<- c(1,12,2,26,3,9,4,23,5,5,6,19)
reclasssoilmat<- matrix(reclasssoil, ncol = 2, byrow = TRUE)
reclasssoilmat## [,1] [,2]
## [1,] 1 12
## [2,] 2 26
## [3,] 3 9
## [4,] 4 23
## [5,] 5 5
## [6,] 6 19reclassreproj<- reclassify(masked60, reclasssoilmat)
plot(reclassreproj, main = "SSURGO DC Texture Dual Depth Method", zlim = c(1,24), col = viridis(24))
writeRaster(reclassreproj, file = 'grassexport/cwt_ws32/ssurgosoil.tif', overwrite = TRUE, options=c("COMPRESS=LZW"))
#plot(SoilsMap2, main = "SSURGO USACE Texture Method")
##this has been modified for quick implementation of shallow-deep classes
##original before depth class subsititution
#reclassdsm<- c(1,9,2,3,3,12,4,0,5,0,6,0,7,0,8,0,9,0,10,0,11,0,12,0,13,0,14,0)
#new
reclassdsm<- c(1,12,2,26,3,9,4,23,5,5,6,19)
reclassdsmmat<- matrix(reclassdsm, ncol = 2, byrow = TRUE)
reclassdsmmat## [,1] [,2]
## [1,] 1 12
## [2,] 2 26
## [3,] 3 9
## [4,] 4 23
## [5,] 5 5
## [6,] 6 19reclassmaskeddsm <- reclassify(maskeddsm, reclassdsmmat)
plot(reclassmaskeddsm, main = "DSM DC Texture Dual Depth Method", zlim = c(1,24), col = viridis(24))
writeRaster(reclassmaskeddsm, file = 'grassexport/cwt_ws32/dsmsoil.tif', overwrite = TRUE, options=c("COMPRESS=LZW"))
#plot(maskeddsm)
reclassstatic<- c(1,12,2,12,3,12,4,12,5,12,6,12,7,12,8,12,9,12,10,12,11,12,12,12,13,12,14,12)
reclassstaticmat<- matrix(reclassstatic, ncol = 2, byrow = TRUE)
reclassstaticmat## [,1] [,2]
## [1,] 1 12
## [2,] 2 12
## [3,] 3 12
## [4,] 4 12
## [5,] 5 12
## [6,] 6 12
## [7,] 7 12
## [8,] 8 12
## [9,] 9 12
## [10,] 10 12
## [11,] 11 12
## [12,] 12 12
## [13,] 13 12
## [14,] 14 12reclassmaskedstatic <- reclassify(maskeddsm, reclassstaticmat)
plot(reclassmaskedstatic, main = "Static Loam Texture", zlim = c(1,24), col = viridis(24))
Combined Soil Plots
library(cowplot)
rssraster<- reclassmaskeddsm
ssurgoraster<- reclassreproj
staticraster<- reclassmaskedstatic
rssrasterdf<- as.data.frame(rssraster, xy = TRUE)
ssurgorasterdf<- as.data.frame(ssurgoraster, xy = TRUE)
staticrasterdf<- as.data.frame(staticraster, xy = TRUE)
rssrasterdf$rss.soiltype.class<- factor(rssrasterdf$rss.soiltype.class)
ssurgorasterdf$ssurgo.soiltype.class<- factor(ssurgorasterdf$ssurgo.soiltype.class)
staticrasterdf$rss.soiltype.class<- factor(staticrasterdf$rss.soiltype.class)
plot3<- ggplot()+
geom_raster(data = rssrasterdf, aes(x = x, y = y, fill = rss.soiltype.class))+
theme_minimal()
plot2<- ggplot()+
geom_raster(data = ssurgorasterdf, aes(x = x, y = y, fill = ssurgo.soiltype.class))+
theme_minimal()
plot1<- ggplot()+
geom_raster(data = staticrasterdf, aes(x = x, y = y, fill = rss.soiltype.class))+
theme_minimal()
combined_soilplot<- plot_grid(plot1, plot2, plot3, ncol =3)
plot(combined_soilplot)
Download Def Files
setwd(system.file("extdata/", package = "RHESSysIOinR"))
#gert::git_clone(url = "https://github.com/RHESSys/ParameterLibrary", path = "defs/rhessys", branch = "master")
listofdefs<- list.files("defs/rhessys/", ".def$", recursive = "TRUE")
file.copy(from = paste0("defs/rhessys/",c(listofdefs)), to = "defs/", overwrite = "TRUE")## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUENLCD Land Cover Classification Legend
11 Open Water 12 Perennial Ice/ Snow 21 Developed, Open Space 22 Developed, Low Intensity 23 Developed, Medium Intensity 24 Developed, High Intensity 31 Barren Land (Rock/Sand/Clay) 41 Deciduous Forest 42 Evergreen Forest 43 Mixed Forest 51 Dwarf Scrub 52 Shrub/Scrub 71 Grassland/Herbaceous 72 Sedge/Herbaceous 73 Lichens 74 Moss 81 Pasture/Hay 82 Cultivated Crops 90 Woody Wetlands 95 Emergent Herbaceous Wetlands
- Alaska only
Because default RHESSys definition files are being used we will reclassify some of this land cover classification for our landscape
41 Deciduous Forest <- 3 veg_deciduous 42 Evergreen Forest <- 4 veg_evergreen 43 Mixed Forest <- 4 veg_evergreen
setwd(system.file("extdata/", package = "RHESSysIOinR"))
LANDCOVER<- get_nlcd(template = x.buff.gcs, label = 'coweeta', dataset = "landcover")
#CANOPY<- get_nlcd(template = x.buff.gcs, label = 'coweeta', year = 2016, dataset = "canopy")
IMPERVIOUS<- get_nlcd(template = x.buff.gcs, label = 'coweeta', dataset = "impervious")
plot(LANDCOVER, main = "NLCD Landcover")
#plot(CANOPY)
plot(IMPERVIOUS, main = "NLCD Impervious")
LANDCOVERREPROJ <- projectRaster(from = LANDCOVER, crs = proj4string(ws32polygon), method = 'ngb')
#CANOPYREPROJ <- projectRaster(from = CANOPY, crs = proj4string(ws32polygon), method = 'ngb')
IMPERVIOUSREPROJ <- projectRaster(from = IMPERVIOUS, crs = proj4string(ws32polygon), method = 'ngb')
## prepare 10m resolution raster and resample prepared soil data to 10m
resample10<- raster(resolution=c(10,10), crs = proj4string(ws32polygon), ext=extent(ws32polygon))
LANDCOVERREPROJRES <- resample(LANDCOVERREPROJ,resample10, method = 'ngb')
{plot(LANDCOVERREPROJRES, main = "Reprojected and Resampled Texture Raster with modified Depth Classes and Watersheds")
plot(st_geometry(x), add = TRUE, col = NA)}
landcovercrop <- crop(LANDCOVERREPROJRES,ws32polygon)
#canopycrop <- crop(CANOPYREPROJ,ws32polygon)
imperviouscrop <- crop(IMPERVIOUSREPROJ,ws32polygon)
maskedlandcover <- mask(landcovercrop,ws32polygon)
#maskedcanopy <- mask(canopycrop,ws32polygon)
maskedimpervious <- mask(imperviouscrop,ws32polygon)
plot(maskedlandcover, main = "Masked Landcover")
#plot(maskedcanopy)
plot(maskedimpervious, main = "Masked Impervious")
{plot(maskedlandcover, main = "Masked Landcover with WS Polygon")
plot(ws32polygon, add = TRUE, col = NA)}
{plot(maskedimpervious, main = "Masked Impervious with WS Polygon")
plot(ws32polygon, add = TRUE, col = NA, border = "white")}
extent(LANDCOVER)## class : Extent
## xmin : 1128045
## xmax : 1136025
## ymin : 1402515
## ymax : 1411275extent(LANDCOVERREPROJ)## class : Extent
## xmin : 271758.1
## xmax : 281385.9
## ymin : 3875918
## ymax : 3886032extent(maskedlandcover)## class : Extent
## xmin : 275066.7
## xmax : 276116.7
## ymin : 3880940
## ymax : 3881570extent(ws32polygon)## class : Extent
## xmin : 275066.7
## xmax : 276116.7
## ymin : 3880940
## ymax : 3881570reclassveg<- c(11,3,12,3,21,3,22,3,23,3,24,3,31,3,41,3,42,4,43,4,51,3,52,3,71,3,72,3,73,3,74,3,81,3,82,3,90,3,95,3)
reclassvegmat<- matrix(reclassveg, ncol = 2, byrow = TRUE)
reclassvegmat## [,1] [,2]
## [1,] 11 3
## [2,] 12 3
## [3,] 21 3
## [4,] 22 3
## [5,] 23 3
## [6,] 24 3
## [7,] 31 3
## [8,] 41 3
## [9,] 42 4
## [10,] 43 4
## [11,] 51 3
## [12,] 52 3
## [13,] 71 3
## [14,] 72 3
## [15,] 73 3
## [16,] 74 3
## [17,] 81 3
## [18,] 82 3
## [19,] 90 3
## [20,] 95 3reclassmaskedveg <- reclassify(maskedlandcover, reclassvegmat)
plot(reclassmaskedveg, main = "Reclassified Vegetation", zlim = c(1,14), col = viridis(14))
writeRaster(reclassmaskedveg, file = 'grassexport/cwt_ws32/nlcdveg.tif', overwrite = TRUE, options=c("COMPRESS=LZW"))Prescribe vegetation rooting depth
setwd(system.file("extdata/", package = "RHESSysIOinR"))
##max root depth seems to be in meters, not clearly defined but inferred from 2021 publication
maxrootdepth = 1
maxrootdepth<- as.data.frame(maxrootdepth)
colnames(maxrootdepth)<-"epc.max_root_depth"
change_def_file("defs/veg_deciduous.def", maxrootdepth, file_name_ext = NULL)## [1] "defs/veg_deciduous/veg_deciduous.def"change_def_file("defs/veg_evergreen.def", maxrootdepth, file_name_ext = NULL)## [1] "defs/veg_evergreen/veg_evergreen.def"LAI data from MODIS
modified from https://cran.r-project.org/web/packages/MODISTools/vignettes/modistools-vignette.html
setwd(system.file("extdata/", package = "RHESSysIOinR"))
library(MODISTools)
bands <- mt_bands(product = "MOD15A2H")
head(bands)## band description units
## 1 FparExtra_QC Extra detail Quality for LAI and FPAR class-flag
## 2 FparLai_QC Quality for LAI and FPAR class-flag
## 3 FparStdDev_500m Standard deviation for FPAR percent
## 4 Fpar_500m Fraction of photosynthetically active radiation percent
## 5 LaiStdDev_500m Standard deviation for LAI m^2/m^2
## 6 Lai_500m Leaf area index m^2/m^2
## valid_range fill_value scale_factor add_offset
## 1 0 to 254 255 <NA> <NA>
## 2 0 to 254 255 <NA> <NA>
## 3 0 to 100 255 0.01 0
## 4 0 to 100 255 0.01 0
## 5 0 to 100 255 0.1 0
## 6 0 to 100 255 0.1 0dates <- mt_dates(product = "MOD15A2H", lat = 35 , lon = -83 )
paged_table(dates)coweeta_lai <- mt_subset (product = "MOD15A2H", lat = 35.06, lon = -83.43,
band = "Lai_500m", start = "2017-01-01", end = "2018-06-01",
km_lr = 20, km_ab = 20, site_name="coweeta",
internal = TRUE, progress = FALSE)
coweeta_lc <- mt_subset (product = "MCD12Q1", lat = 35.06, lon = -83.43,
band = "LC_Type1", start = "2017-01-01", end = "2018-06-01",
km_lr = 20, km_ab = 20, site_name="coweeta",
internal = TRUE, progress = FALSE)
head(coweeta_lai)## xllcorner yllcorner cellsize nrows ncols band units scale
## 1.1 -7612691.23 3879780.73 463.312716528 81 81 Lai_500m m^2/m^2 0.1
## 2.1 -7612691.23 3879780.73 463.312716528 81 81 Lai_500m m^2/m^2 0.1
## 3.1 -7612691.23 3879780.73 463.312716528 81 81 Lai_500m m^2/m^2 0.1
## 4.1 -7612691.23 3879780.73 463.312716528 81 81 Lai_500m m^2/m^2 0.1
## 5.1 -7612691.23 3879780.73 463.312716528 81 81 Lai_500m m^2/m^2 0.1
## 6.1 -7612691.23 3879780.73 463.312716528 81 81 Lai_500m m^2/m^2 0.1
## latitude longitude site product start end complete
## 1.1 35.06 -83.43 coweeta MOD15A2H 2017-01-01 2018-06-01 TRUE
## 2.1 35.06 -83.43 coweeta MOD15A2H 2017-01-01 2018-06-01 TRUE
## 3.1 35.06 -83.43 coweeta MOD15A2H 2017-01-01 2018-06-01 TRUE
## 4.1 35.06 -83.43 coweeta MOD15A2H 2017-01-01 2018-06-01 TRUE
## 5.1 35.06 -83.43 coweeta MOD15A2H 2017-01-01 2018-06-01 TRUE
## 6.1 35.06 -83.43 coweeta MOD15A2H 2017-01-01 2018-06-01 TRUE
## modis_date calendar_date tile proc_date pixel value
## 1.1 A2017001 2017-01-01 h11v05 2021252144002 1 5
## 2.1 A2017009 2017-01-09 h11v05 2021257201833 1 3
## 3.1 A2017017 2017-01-17 h11v05 2021258170909 1 7
## 4.1 A2017025 2017-01-25 h11v05 2021264192656 1 8
## 5.1 A2017033 2017-02-02 h11v05 2021265175312 1 4
## 6.1 A2017041 2017-02-10 h11v05 2021266102506 1 7coweetamodis <- coweeta_lc %>%
rename("lc" = "value") %>%
select("lc","pixel") %>%
right_join(coweeta_lai, by = "pixel")
coweetamodis <- coweetamodis %>%
filter(value <= 100,
lc %in% c("1","5")) %>%
mutate(lc = ifelse(lc == 1, "ENF","DBF")) %>%
group_by(lc, calendar_date) %>%
summarize(doy = as.numeric(format(as.Date(calendar_date)[1],"%j")),
lai_mean = median(value * as.double(scale)))
ggplot(coweetamodis, aes(x = doy, y = lai_mean)) +
geom_point() +
geom_smooth(span = 0.3, method = "loess") +
labs(x = "day of year (DOY)",
y = "leaf area index (LAI)") +
theme_minimal() +
facet_wrap(~ lc)
coweeta_lai <- mt_subset (product = "MOD15A2H", lat = 35.06, lon = -83.43,
band = "Lai_500m", start = "2017-05-20", end = "2017-06-01",
km_lr = 20, km_ab = 20, site_name="coweeta",
internal = TRUE, progress = FALSE)
LAI_raster <- mt_to_raster(df=coweeta_lai, reproject = TRUE)
x <- read_sf(watershedshapefile)
x.buff <- st_as_sf(st_buffer(st_as_sfc(st_bbox(x)), dist = 1000))
x.gcs <- st_transform(x, 4326)
x.buff.gcs <- st_transform(x.buff, 4326)
{plot(LAI_raster)
plot(st_geometry(x.buff.gcs), add = TRUE, col = NA, lwd = 2)
plot(st_geometry(x.gcs), add = TRUE, col = NA)
}
Create Templates for RHESSYs
setwd(system.file("extdata/", package = "RHESSysIOinR"))
# template_read(template)
# would be useful to modify or read RHESSys template in R, currently template_read does not display nicely
# Create Templates for each Treatment
dir.create("templates")
{templatetext<- "1
../defs/basin.def
1
../defs/hillslope.def
1
../defs/zone.def
14
../defs/soil_clay.def
../defs/soil_clayloam.def
../defs/soil_loam.def
../defs/soil_loamysand.def
../defs/soil_rock.def
../defs/soil_sand.def
../defs/soil_sandyclay.def
../defs/soil_sandyclayloam.def
../defs/soil_sandyloam.def
../defs/soil_silt.def
../defs/soil_siltyclay.def
../defs/soil_siltyclayloam.def
../defs/soil_siltyloam.def
../defs/soil_water.def
../defs/soil_shallowloam.def
../defs/soil_shallowsandyclayloam
../defs/soil_shallowsandyloam
1
../defs/lu_undev.def
9
../defs/veg_deciduous/veg_deciduous.def
../defs/veg_evergreen/veg_evergreen.def
../defs/veg_deciduous_BES.def
../defs/veg_eucalypt.def
../defs/veg_grass.def
../defs/veg_lawn_2cm.def
../defs/veg_lawn_5cm.def
../defs/veg_lawn_10cm.def
../defs/veg_nonveg.def
1
../clim/cwtws32extended.base
_world basin_ws32.tif 1
_basin basin_ws32.tif 1
x value 0.0
y value 0.0
z aver dem10.tif
basin_parm_ID dvalue 1
latitude value 35.04
basin_n_basestations dvalue 0
_hillslope subbasins.tif 1
x value 0.0
y value 0.0
z aver dem10.tif
hill_parm_ID dvalue 1
gw.storage value 0.0
gw.NO3 value 0.0
hillslope_n_basestations dvalue 0
_zone patch.tif 1
x value 0.0
y value 0.0
z aver dem10.tif
zone_parm_ID dvalue 1
area area
slope aver slope10.tif
aspect spavg aspect10.tif slope10.tif
precip_lapse_rate value 1.0
e_horizon eqn 0.001 0 e_horizon_1000.tif
w_horizon eqn 0.001 0 w_horizon_1000.tif
zone_n_basestations dvalue 1
zone_basestation_ID dvalue 101
_patch patch.tif 1
x value 0.0
y value 0.0
z aver dem10.tif
soil_parm_ID mode ssurgosoil.tif
landuse_parm_ID dvalue 1
fire_parm_ID value 1.0
area area
slope aver slope10.tif
lna eqn 0.001 0 topidx10_100.tif
Ksat_vertical value 1.0
mpar value 0.12
rz_storage value 0.0
unsat_storage value 0.0
sat_deficit value 0.0
snowpack.water_equivalent_depth value 0.28
snowpack.water_depth value 0.0
snowpack.T value -10.0
snowpack.surface_age value 0.0
snowpack.energy_deficit value -0.5
litter.cover_fraction value 1.0
litter.rain_stored value 0.00001544
litter_cs.litr1c value 0.00158273
litter_ns.litr1n value 0.00007521
litter_cs.litr2c value 0.01056238
litter_cs.litr3c value 0.00094835
litter_cs.litr4c value 0.00641051
soil_cs.soil1c value 0.00428524
soil_ns.sminn value 0.00002811
soil_ns.nitrate value 0.00003044
soil_cs.soil2c value 0.00304513
soil_cs.soil3c value 0.05220800
soil_cs.soil4c value 0.37839627
patch_n_basestations dvalue 0
_canopy_strata patch.tif 1
veg_parm_ID mode nlcdveg.tif
cover_fraction value 1.0
gap_fraction value 0.0
rootzone.depth value 1.0
snow_stored value 0.0
rain_stored value 0.0
cs.pool value 0.0
cs.leafc value 0.0
cs.dead_leafc value 0.0
cs.leafc_store value 0.0
cs.leafc_transfer value 0.0
cs.live_stemc value 0.0
cs.livestemc_store value 0.0
cs.live.stemc_transfer value 0.0
cs.dead_stem value 0.0
cs.deadstemc_store value 0.0
cs.deadstemc_transfer value 0.0
cs.live_crootc value 0.0
cs.livecrootc_store value 0.0
cs.livecrootc_transfer value 0.0
cs.dead_crootc value 2.0
cs.deadcrootc_store value 0.0
cs.deadcrootc_transfer value 0.0
cs.frootc value 0.0
cs.frootc_store value 0.0
cs.frootc_transfer value 0.0
cs_cwdc value 0.0
epv.prev_leafcalloc value 0.0
ns.pool value 0.0
ns.leafn value 0.0
ns.dead_leafn value 0.0
ns.leafn_store value 0.0
ns.leafn_transfer value 0.0
ns.live_stemn value 0.0
ns.livestemn_store value 0.0
ns.livestemn_transfer value 0.0
ns.dead_stem value 0.0
ns.deadstemn_store value 0.0
ns.deadstemn_transfer value 0.0
ns.live_crootn value 0.0
ns.livecrootn_store value 0.0
ns.livecrootn_transfer value 0.0
ns.dead_crootn value 0.0
ns.deadcrootn_store value 0.0
ns.deadcrootn_transfer value 0.0
ns.frootn value 0.0
ns.frootn_store value 0.0
ns.frootn_transfer value 0.0
ns.cwdn value 0.0
ns.retransn value 0.0
epv.wstress_days dvalue 0
epv.min_fparabs value 0.0
epv.min_vwc value 0.0
canopy_strata_n_basestations dvalue 0
"}
write(templatetext,file = "templates/RHESSYSinR_Coweeta_WS32_Spinup.template")
#Static
{templatestatic<- "1
../defs/basin.def
1
../defs/hillslope.def
1
../defs/zone.def
14
../defs/soil_clay.def
../defs/soil_clayloam.def
../defs/soil_loam.def
../defs/soil_loamysand.def
../defs/soil_rock.def
../defs/soil_sand.def
../defs/soil_sandyclay.def
../defs/soil_sandyclayloam.def
../defs/soil_sandyloam.def
../defs/soil_silt.def
../defs/soil_siltyclay.def
../defs/soil_siltyclayloam.def
../defs/soil_siltyloam.def
../defs/soil_water.def
../defs/soil_shallowloam.def
../defs/soil_shallowsandyclayloam
../defs/soil_shallowsandyloam
1
../defs/lu_undev.def
2
../defs/veg_deciduous/veg_deciduous.def
../defs/veg_evergreen/veg_evergreen.def
1
../clim/cwtws32extended.base
_world basin_ws32.tif 1
_basin basin_ws32.tif 1
x value 0.0
y value 0.0
z aver dem10.tif
basin_parm_ID dvalue 1
latitude value 35.04
basin_n_basestations dvalue 0
_hillslope subbasins.tif 1
x value 0.0
y value 0.0
z aver dem10.tif
hill_parm_ID dvalue 1
gw.storage value 0.0
gw.NO3 value 0.0
hillslope_n_basestations dvalue 0
_zone patch.tif 1
x value 0.0
y value 0.0
z aver dem10.tif
zone_parm_ID dvalue 1
area area
slope aver slope10.tif
aspect spavg aspect10.tif slope10.tif
precip_lapse_rate value 1.0
e_horizon eqn 0.001 0 e_horizon_1000.tif
w_horizon eqn 0.001 0 w_horizon_1000.tif
zone_n_basestations dvalue 1
zone_basestation_ID dvalue 101
_patch patch.tif 1
x value 0.0
y value 0.0
z aver dem10.tif
soil_parm_ID dvalue 9
landuse_parm_ID dvalue 1
fire_parm_ID value 1.0
area area
slope aver slope10.tif
lna eqn 0.001 0 topidx10_100.tif
Ksat_vertical value 1.0
mpar value 0.12
rz_storage value 0.0
unsat_storage value 0.0
sat_deficit value 0.0
snowpack.water_equivalent_depth value 0.28
snowpack.water_depth value 0.0
snowpack.T value -10.0
snowpack.surface_age value 0.0
snowpack.energy_deficit value -0.5
litter.cover_fraction value 1.0
litter.rain_stored value 0.00001544
litter_cs.litr1c value 0.00158273
litter_ns.litr1n value 0.00007521
litter_cs.litr2c value 0.01056238
litter_cs.litr3c value 0.00094835
litter_cs.litr4c value 0.00641051
soil_cs.soil1c value 0.00428524
soil_ns.sminn value 0.00002811
soil_ns.nitrate value 0.00003044
soil_cs.soil2c value 0.00304513
soil_cs.soil3c value 0.05220800
soil_cs.soil4c value 0.37839627
patch_n_basestations dvalue 0
_canopy_strata patch.tif 1
veg_parm_ID mode nlcdveg.tif
cover_fraction value 1.0
gap_fraction value 0.0
rootzone.depth value 1.0
snow_stored value 0.0
rain_stored value 0.0
cs.pool value 0.0
cs.leafc value 0.0
cs.dead_leafc value 0.0
cs.leafc_store value 0.0
cs.leafc_transfer value 0.0
cs.live_stemc value 0.0
cs.livestemc_store value 0.0
cs.live.stemc_transfer value 0.0
cs.dead_stem value 0.0
cs.deadstemc_store value 0.0
cs.deadstemc_transfer value 0.0
cs.live_crootc value 0.0
cs.livecrootc_store value 0.0
cs.livecrootc_transfer value 0.0
cs.dead_crootc value 2.0
cs.deadcrootc_store value 0.0
cs.deadcrootc_transfer value 0.0
cs.frootc value 0.0
cs.frootc_store value 0.0
cs.frootc_transfer value 0.0
cs_cwdc value 0.0
epv.prev_leafcalloc value 0.0
ns.pool value 0.0
ns.leafn value 0.0
ns.dead_leafn value 0.0
ns.leafn_store value 0.0
ns.leafn_transfer value 0.0
ns.live_stemn value 0.0
ns.livestemn_store value 0.0
ns.livestemn_transfer value 0.0
ns.dead_stem value 0.0
ns.deadstemn_store value 0.0
ns.deadstemn_transfer value 0.0
ns.live_crootn value 0.0
ns.livecrootn_store value 0.0
ns.livecrootn_transfer value 0.0
ns.dead_crootn value 0.0
ns.deadcrootn_store value 0.0
ns.deadcrootn_transfer value 0.0
ns.frootn value 0.0
ns.frootn_store value 0.0
ns.frootn_transfer value 0.0
ns.cwdn value 0.0
ns.retransn value 0.0
epv.wstress_days dvalue 0
epv.min_fparabs value 0.0
epv.min_vwc value 0.0
canopy_strata_n_basestations dvalue 0
"}
#SSURGO
{templatessurgo<- "1
../defs/basin.def
1
../defs/hillslope.def
1
../defs/zone.def
14
../defs/soil_clay.def
../defs/soil_clayloam.def
../defs/soil_loam.def
../defs/soil_loamysand.def
../defs/soil_rock.def
../defs/soil_sand.def
../defs/soil_sandyclay.def
../defs/soil_sandyclayloam.def
../defs/soil_sandyloam.def
../defs/soil_silt.def
../defs/soil_siltyclay.def
../defs/soil_siltyclayloam.def
../defs/soil_siltyloam.def
../defs/soil_water.def
../defs/soil_shallowloam.def
../defs/soil_shallowsandyclayloam
../defs/soil_shallowsandyloam
1
../defs/lu_undev.def
2
../defs/veg_deciduous/veg_deciduous.def
../defs/veg_evergreen/veg_evergreen.def
1
../clim/cwtws32extended.base
_world basin_ws32.tif 1
_basin basin_ws32.tif 1
x value 0.0
y value 0.0
z aver dem10.tif
basin_parm_ID dvalue 1
latitude value 35.04
basin_n_basestations dvalue 0
_hillslope subbasins.tif 1
x value 0.0
y value 0.0
z aver dem10.tif
hill_parm_ID dvalue 1
gw.storage value 0.0
gw.NO3 value 0.0
hillslope_n_basestations dvalue 0
_zone patch.tif 1
x value 0.0
y value 0.0
z aver dem10.tif
zone_parm_ID dvalue 1
area area
slope aver slope10.tif
aspect spavg aspect10.tif slope10.tif
precip_lapse_rate value 1.0
e_horizon eqn 0.001 0 e_horizon_1000.tif
w_horizon eqn 0.001 0 w_horizon_1000.tif
zone_n_basestations dvalue 1
zone_basestation_ID dvalue 101
_patch patch.tif 1
x value 0.0
y value 0.0
z aver dem10.tif
soil_parm_ID mode ssurgosoil.tif
landuse_parm_ID dvalue 1
fire_parm_ID value 1.0
area area
slope aver slope10.tif
lna eqn 0.001 0 topidx10_100.tif
Ksat_vertical value 1.0
mpar value 0.12
rz_storage value 0.0
unsat_storage value 0.0
sat_deficit value 0.0
snowpack.water_equivalent_depth value 0.28
snowpack.water_depth value 0.0
snowpack.T value -10.0
snowpack.surface_age value 0.0
snowpack.energy_deficit value -0.5
litter.cover_fraction value 1.0
litter.rain_stored value 0.00001544
litter_cs.litr1c value 0.00158273
litter_ns.litr1n value 0.00007521
litter_cs.litr2c value 0.01056238
litter_cs.litr3c value 0.00094835
litter_cs.litr4c value 0.00641051
soil_cs.soil1c value 0.00428524
soil_ns.sminn value 0.00002811
soil_ns.nitrate value 0.00003044
soil_cs.soil2c value 0.00304513
soil_cs.soil3c value 0.05220800
soil_cs.soil4c value 0.37839627
patch_n_basestations dvalue 0
_canopy_strata patch.tif 1
veg_parm_ID mode nlcdveg.tif
cover_fraction value 1.0
gap_fraction value 0.0
rootzone.depth value 1.0
snow_stored value 0.0
rain_stored value 0.0
cs.pool value 0.0
cs.leafc value 0.0
cs.dead_leafc value 0.0
cs.leafc_store value 0.0
cs.leafc_transfer value 0.0
cs.live_stemc value 0.0
cs.livestemc_store value 0.0
cs.live.stemc_transfer value 0.0
cs.dead_stem value 0.0
cs.deadstemc_store value 0.0
cs.deadstemc_transfer value 0.0
cs.live_crootc value 0.0
cs.livecrootc_store value 0.0
cs.livecrootc_transfer value 0.0
cs.dead_crootc value 2.0
cs.deadcrootc_store value 0.0
cs.deadcrootc_transfer value 0.0
cs.frootc value 0.0
cs.frootc_store value 0.0
cs.frootc_transfer value 0.0
cs_cwdc value 0.0
epv.prev_leafcalloc value 0.0
ns.pool value 0.0
ns.leafn value 0.0
ns.dead_leafn value 0.0
ns.leafn_store value 0.0
ns.leafn_transfer value 0.0
ns.live_stemn value 0.0
ns.livestemn_store value 0.0
ns.livestemn_transfer value 0.0
ns.dead_stem value 0.0
ns.deadstemn_store value 0.0
ns.deadstemn_transfer value 0.0
ns.live_crootn value 0.0
ns.livecrootn_store value 0.0
ns.livecrootn_transfer value 0.0
ns.dead_crootn value 0.0
ns.deadcrootn_store value 0.0
ns.deadcrootn_transfer value 0.0
ns.frootn value 0.0
ns.frootn_store value 0.0
ns.frootn_transfer value 0.0
ns.cwdn value 0.0
ns.retransn value 0.0
epv.wstress_days dvalue 0
epv.min_fparabs value 0.0
epv.min_vwc value 0.0
canopy_strata_n_basestations dvalue 0
"}
#DSM
{templatedsm<- "1
../defs/basin.def
1
../defs/hillslope.def
1
../defs/zone.def
14
../defs/soil_clay.def
../defs/soil_clayloam.def
../defs/soil_loam.def
../defs/soil_loamysand.def
../defs/soil_rock.def
../defs/soil_sand.def
../defs/soil_sandyclay.def
../defs/soil_sandyclayloam.def
../defs/soil_sandyloam.def
../defs/soil_silt.def
../defs/soil_siltyclay.def
../defs/soil_siltyclayloam.def
../defs/soil_siltyloam.def
../defs/soil_water.def
../defs/soil_shallowloam.def
../defs/soil_shallowsandyclayloam
../defs/soil_shallowsandyloam
1
../defs/lu_undev.def
1
../defs/veg_deciduous/veg_deciduous.def
../defs/veg_evergreen/veg_evergreen.def
1
../clim/cwtws32extended.base
_world basin_ws32.tif 1
_basin basin_ws32.tif 1
x value 0.0
y value 0.0
z aver dem10.tif
basin_parm_ID dvalue 1
latitude value 35.04
basin_n_basestations dvalue 0
_hillslope subbasins.tif 1
x value 0.0
y value 0.0
z aver dem10.tif
hill_parm_ID dvalue 1
gw.storage value 0.0
gw.NO3 value 0.0
hillslope_n_basestations dvalue 0
_zone patch.tif 1
x value 0.0
y value 0.0
z aver dem10.tif
zone_parm_ID dvalue 1
area area
slope aver slope10.tif
aspect spavg aspect10.tif slope10.tif
precip_lapse_rate value 1.0
e_horizon eqn 0.001 0 e_horizon_1000.tif
w_horizon eqn 0.001 0 w_horizon_1000.tif
zone_n_basestations dvalue 1
zone_basestation_ID dvalue 101
_patch patch.tif 1
x value 0.0
y value 0.0
z aver dem10.tif
soil_parm_ID mode dsmsoil.tif
landuse_parm_ID dvalue 1
fire_parm_ID value 1.0
area area
slope aver slope10.tif
lna eqn 0.001 0 topidx10_100.tif
Ksat_vertical value 1.0
mpar value 0.12
rz_storage value 0.0
unsat_storage value 0.0
sat_deficit value 0.0
snowpack.water_equivalent_depth value 0.28
snowpack.water_depth value 0.0
snowpack.T value -10.0
snowpack.surface_age value 0.0
snowpack.energy_deficit value -0.5
litter.cover_fraction value 1.0
litter.rain_stored value 0.00001544
litter_cs.litr1c value 0.00158273
litter_ns.litr1n value 0.00007521
litter_cs.litr2c value 0.01056238
litter_cs.litr3c value 0.00094835
litter_cs.litr4c value 0.00641051
soil_cs.soil1c value 0.00428524
soil_ns.sminn value 0.00002811
soil_ns.nitrate value 0.00003044
soil_cs.soil2c value 0.00304513
soil_cs.soil3c value 0.05220800
soil_cs.soil4c value 0.37839627
patch_n_basestations dvalue 0
_canopy_strata patch.tif 1
veg_parm_ID mode nlcdveg.tif
cover_fraction value 1.0
gap_fraction value 0.0
rootzone.depth value 1.0
snow_stored value 0.0
rain_stored value 0.0
cs.pool value 0.0
cs.leafc value 0.0
cs.dead_leafc value 0.0
cs.leafc_store value 0.0
cs.leafc_transfer value 0.0
cs.live_stemc value 0.0
cs.livestemc_store value 0.0
cs.live.stemc_transfer value 0.0
cs.dead_stem value 0.0
cs.deadstemc_store value 0.0
cs.deadstemc_transfer value 0.0
cs.live_crootc value 0.0
cs.livecrootc_store value 0.0
cs.livecrootc_transfer value 0.0
cs.dead_crootc value 2.0
cs.deadcrootc_store value 0.0
cs.deadcrootc_transfer value 0.0
cs.frootc value 0.0
cs.frootc_store value 0.0
cs.frootc_transfer value 0.0
cs_cwdc value 0.0
epv.prev_leafcalloc value 0.0
ns.pool value 0.0
ns.leafn value 0.0
ns.dead_leafn value 0.0
ns.leafn_store value 0.0
ns.leafn_transfer value 0.0
ns.live_stemn value 0.0
ns.livestemn_store value 0.0
ns.livestemn_transfer value 0.0
ns.dead_stem value 0.0
ns.deadstemn_store value 0.0
ns.deadstemn_transfer value 0.0
ns.live_crootn value 0.0
ns.livecrootn_store value 0.0
ns.livecrootn_transfer value 0.0
ns.dead_crootn value 0.0
ns.deadcrootn_store value 0.0
ns.deadcrootn_transfer value 0.0
ns.frootn value 0.0
ns.frootn_store value 0.0
ns.frootn_transfer value 0.0
ns.cwdn value 0.0
ns.retransn value 0.0
epv.wstress_days dvalue 0
epv.min_fparabs value 0.0
epv.min_vwc value 0.0
canopy_strata_n_basestations dvalue 0
"}
write(templatestatic,file = "templates/RHESSYSinR_Coweeta_WS32_static.template")
write(templatetext,file = "templates/RHESSYSinR_Coweeta_WS32_ssurgo.template")
write(templatedsm,file = "templates/RHESSYSinR_Coweeta_WS32_dsm.template")Create Tec File
setwd(system.file("extdata/", package = "RHESSysIOinR"))
{tecfiletext<-
"2015 11 1 1 print_daily_on
2015 11 1 2 print_daily_growth_on
"
write(tecfiletext,file = "tecfiles/tec_daily")}Locate RHESSys Spatial Inputs and Specify Template to be used
Create RHESSys Worldfile and Flowtable for 7.4 based on Spatial Inputs and Template
Raster files and template information are used to setup RHESSys worldfile. Flowtable created by “RHESSysPreprocess” must be converted to work. Unsure why. Preprocessing is controlled by template file.
Template file ‘CWWS32_STATIC.template’ is used to run model normally template file ‘CWWS32_STATICspinup.template’ is used to spin up model. Only difference is artificially extended climate dataset.
must setup automatic spinup and preprocessing for all treatments, right now spinup requires manual template input and 2 runs for creation
setwd(system.file("extdata/", package = "RHESSysIOinR"))
# Optional Flowtable Spatial Data
roads = "roads.tif"
# impervious = "impervious_map"
# roofs = "roofs_map"
# With certain configurations roads does not work. Returns message: road width cannot be 0 during preprocessing
#
# RHESSysPreprocess(
# template = "templates/RHESSYSinR_Coweeta_WS32_dsm.template",
# name = "CWWS32",
# map_dir = "grassexport/cwt_ws32",
# streams = "streams.tif",
# overwrite = TRUE,
# header = TRUE,
# parallel = FALSE,
# make_stream = TRUE)
## Make sure to have parallel turned off when preprocessing, convert flowtables after creation for compatibility
#convert_flowtable("CWWS32.flow","CWWS321.flow",make_stream=TRUE)
#Try alternate method
#
world_gen(template = "templates/RHESSYSinR_Coweeta_WS32_dsm.template",
worldfile = "CWWS32",
map_dir = "grassexport/cwt_ws32",
overwrite = TRUE,
header = TRUE)## Reading in maps...
## Setting NaNs to NA.
## Trimming NAs.
## Masking maps by 'world' layer map.
## Checking for missing data within bounds of world map.
## Finished reading in maps
## [1] Writing worldfile
##
|
| | 0%
|
|======================= | 33%
|
|=============================================== | 67%
|
|======================================================================| 100%
## [1] Created worldfile: ./CWWS32.world
## [1] Created header file: ./CWWS32.hdr## [[1]]
## MapName Map
## [1,] "world" "basin_ws32.tif"
## [2,] "basin" "basin_ws32.tif"
## [3,] "hillslope" "subbasins.tif"
## [4,] "zone" "patch.tif"
## [5,] "patch" "patch.tif"
## [6,] "strata" "patch.tif"
## [7,] "z" "dem10.tif"
## [8,] "z" "dem10.tif"
## [9,] "z" "dem10.tif"
## [10,] "slope" "slope10.tif"
## [11,] "aspect" "aspect10.tif"
## [12,] "e_horizon" "e_horizon_1000.tif"
## [13,] "w_horizon" "w_horizon_1000.tif"
## [14,] "z" "dem10.tif"
## [15,] "soil_parm_ID" "dsmsoil.tif"
## [16,] "slope" "slope10.tif"
## [17,] "lna" "topidx10_100.tif"
## [18,] "veg_parm_ID" "nlcdveg.tif"
## [19,] "cell_length" "10"
## [20,] "streams" "none"
## [21,] "roads" "none"
## [22,] "impervious" "none"
## [23,] "roofs" "none"
##
## [[2]]
## NULLcreate_flownet(template = "templates/RHESSYSinR_Coweeta_WS32_dsm.template",
flownet_name = "CWWS32",
map_dir = "grassexport/cwt_ws32",
streams = "streams.tif",
overwrite = TRUE,
parallel = FALSE,
make_stream = TRUE)## Reading in maps...
## Setting NaNs to NA.
## Trimming NAs.
## Masking maps by 'world' layer map.
## Checking for missing data within bounds of world map.
## Finished reading in maps
##
## Building flowtable
## Generating unique patch IDs
## Finding patch neighbors
##
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========= | 14%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|============== | 21%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================ | 24%
|
|================= | 24%
|
|================= | 25%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|===================== | 31%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 32%
|
|======================= | 33%
|
|======================= | 34%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|============================== | 44%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 46%
|
|================================= | 47%
|
|================================= | 48%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 49%
|
|=================================== | 50%
|
|=================================== | 51%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 53%
|
|===================================== | 54%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|========================================== | 61%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================ | 64%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|=================================================== | 74%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 79%
|
|======================================================== | 80%
|
|======================================================== | 81%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|========================================================== | 84%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================= | 94%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%
## Buildling flowtable list
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========= | 14%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|============== | 21%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================ | 24%
|
|================= | 24%
|
|================= | 25%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|===================== | 31%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 32%
|
|======================= | 33%
|
|======================= | 34%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|============================== | 44%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 46%
|
|================================= | 47%
|
|================================= | 48%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 49%
|
|=================================== | 50%
|
|=================================== | 51%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 53%
|
|===================================== | 54%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|========================================== | 61%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================ | 64%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|=================================================== | 74%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 79%
|
|======================================================== | 80%
|
|======================================================== | 81%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|========================================================== | 84%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================= | 94%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%
## Writing flowtable
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========= | 14%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|============== | 21%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================ | 24%
|
|================= | 24%
|
|================= | 25%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|===================== | 31%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 32%
|
|======================= | 33%
|
|======================= | 34%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|============================== | 44%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 46%
|
|================================= | 47%
|
|================================= | 48%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 49%
|
|=================================== | 50%
|
|=================================== | 51%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 53%
|
|===================================== | 54%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|========================================== | 61%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================ | 64%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|=================================================== | 74%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 79%
|
|======================================================== | 80%
|
|======================================================== | 81%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|========================================================== | 84%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================= | 94%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%
## Created flowtable: ./CWWS32.flow